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SUMMARY 


A  focused  research  program  has  been  conducted  to  investigate  the  dynamic 
behavior  of  liquid-propellant  droplets  in  supercritical  forced-convective  environments. 
The  purpose  is  to  establish  a  solid  theoretical  basis  for  enhancing  the  understanding  of 
liquid-propellant  droplet  vaporization,  combustion,  and  dynamics  at  supercritical 
conditions,  with  emphasis  placed  on  the  effect  of  forced  convection.  A  variety  of  liquid 
propellants  and  propellant  simulants,  including  hydrocarbon  and  cryogenic  fluids,  at  both 
steady  and  oscillatory  conditions  are  treated  systematically.  The  formulation  is  based  on 
the  full  conservation  equations  for  both  gas  and  liquid  phases,  and  accommodates 
variable  properties  and  finite-rate  chemical  kinetics.  Full  account  is  taken  of 
thermodynamic  non-idealities  and  transport  anomalies  at  high  pressures,  as  well  as 
liquid-vapor  phase  equilibria  for  multi-component  mixtures.  Because  the  model  allows 
for  solutions  from  first  principles,  a  systematic  examination  of  droplet  behavior  over 
wide  ranges  of  temperature  and  pressure  is  made  possible. 

The  work  represents  a  series  of  fundamental  studies  of  liquid-propellant  droplet 
dynamics  and  combustion  at  supercritical  conditions.  Results  have  not  only  enhanced  the 
basic  understanding  of  the  problem,  but  can  also  be  used  to  improve  the  predictive 
capabilities  of  contemporary  rocket-engine  performance  and  stability  codes.  In 
particular,  correlations  of  droplet  mass,  momentum,  and  energy  transfer  rates  over  broad 
ranges  of  thermodynamic  and  flow  conditions  are  established  for  use  in  the  study  of  spray 
combustion.  Dynamic  responses  of  droplet  vaporization  and  combustion  to  ambient 
pressure  and  velocity  oscillations  are  also  treated.  The  resultant  combustion  response 
functions  serve  as  a  primary  factor  in  determining  the  stability  characteristics  of  a  rocket 
engine. 


Throughout  the  program,  close  collaboration  with  various  experimental  groups  at 
the  PSU  Propulsion  Engineering  Research  Center  and  the  Air  Force  Research  Laboratory 
is  maintained,  so  that  optimal  benefit  from  all  efforts  are  obtained. 
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TECHNICAL  DISCUSSION 


Liquid-droplet  vaporization  and  combustion  in  supercritical  environments  have 
long  been  matters  of  serious  practical  concern  in  combustion  technology,  mainly  due  to 
the  necessity  of  developing  high-pressure  combustion  devices  such  as  liquid-propellant 
rocket  engines,  diesel  engines  and  liquid  propellant  guns.  Although  several  studies  have 
been  devoted  to  this  problem,  a  number  of  fundamental  issues  regarding  the  attainability 
of  critical  conditions  and  droplet  gasification  and  burning  mechanisms  at  high  pressures 
remain  unresolved.  Most  of  the  existing  theories  are  based  on  certain  assumptions  and 
empirical  correlations  that  are  extrapolated  from  atmospherical  conditions,  with  their 
validity  for  high-pressure  applications  subject  to  clarification.  Furthermore,  no  unified 
analysis  of  the  entire  droplet  history,  in  particular  the  transition  from  the  subcritcal  to 
supercritical  state,  was  conducted.  Almost  all  of  the  models  for  supercritical  combustion 
have  assumed  a  priori  that  the  droplet  reaches  its  critical  state  instantaneously  upon 
introduction  into  a  supercritical  environment.  Neither  initial  heatup  transients  nor 
nonuniformities  of  liquid-phase  properties  are  considered.  In  addition,  the  treatment  of 
transport  properties  and  thermodynamics  nonidealities  is  overly  simplified  to  yield 
complete  information.  Recent  reviews  of  this  subject  are  given  in  Refs.  1-4. 

The  physical  model  treated  here  is  an  isolated  liquid-propellant  droplet  (or  an 
array  of  droplets)  when  suddenly  confronted  with  a  supercritical  fluid  flow.  The  initial 
temperature  of  the  droplet  is  subcritical.  As  the  droplet  is  heated  by  the  ambient  gas,  its 
temperature  increases  and  finally  exceeds  the  critical  point.  Several  important  aspects 
must  be  noted  during  this  process.  First,  when  the  droplet  surface  approaches  its 
thermodynamic  critical  state,  the  difference  in  density  between  gas  and  liquid  phases 
becomes  smaller.  The  characteristic  times  of  the  transport  processes  near  the  droplet 
surface  in  both  phases  have  the  same  order  of  magnitude.  Therefore,  the  transient  effects 
in  the  gas  phase  are  as  important  as  those  in  the  liquid  phase,  and  the  quasi-steady 
condition  may  never  be  reached  during  the  lifetime  of  the  droplet.  Second,  the  latent  heat 
of  vaporization  decreases  to  zero  at  the  critical  point.  Conventional  low-pressure  models 
may  erroneously  predict  the  vaporization  rate  if  the  variation  of  latent  heat  with  pressure 
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is  not  properly  taken  into  account.  In  addition,  if  the  droplet  is  moving,  the  behavior  of 
liquid  deformation  and  breakup  may  be  altered  considerably  due  to  the  diminished  value 
of  surface  tension.  Third,  at  high  pressures,  effects  of  thermodynamic  non-idealities  and 
property  variations  play  decisive  roles  in  determining  transport  properties  and  interfacial 
thermod3mamic  relationships.  The  solubility  of  the  ambient  gas  in  the  liquid  phase 
increases  with  pressure,  and  the  classical  Raoult’s  law  for  ideal  mixtures  is  not  applicable 
for  phase-equilibrium  analysis.  One  must  develop  a  more  comprehensive  model  for 
vapor-liquid  interface  conditions  in  terms  of  fugacity.  Finally,  when  the  droplet  exceeds 
its  critical  state,  it  essentially  becomes  a  “puff’  of  dense  fluid.  The  entire  field  becomes  a 
continuous  medium,  and  no  distinct  interfacial  boundary  can  be  identified. 

The  primary  purpose  of  this  research  project  is  to  establish  a  theoretical 
framework  within  which  various  issues  of  supercritical  droplet  vaporization  and 
combustion  can  be  addressed  systematically.  Specific  tasks  conducted  include: 

1.  investigation  of  supercritical  droplet  vaporization  and  combustion  in  steady 
convective  environments;  and 

2.  study  of  dynamic  responses  of  droplet  vaporization  and  combustion  to 
ambient  flow  oscillations. 

Major  results  obtained  from  these  tasks  are  summarized  below. 

1.  Supercritical  Droplet  Vaporization  and  Combustion  in  Steady  Convective 
Environments 

Several  comprehensive  theoretical  analyses  have  been  established  to  investigate 
droplet  vaporization  and  combustion  in  forced  convective  environments  at  supercritical 
conditions.  Both  hydrocarbon  and  cryogenic  propellant  droplets  are  considered  over  a 
pressure  range  from  one  to  400  atm.  Specific  processes  treated  in  this  task  are: 

1.  complete  time  history  of  droplet  gasification,  including  the  transition  from  the 
subcritical  to  supercritical  state; 

2.  ignition  and  flame  development; 

3.  droplet  dynamics,  including  deformation,  stripping,  and  shattering;  and 
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4.  interphase  transport  between  droplet  and  surrounding  flow. 


Figure  1  shows  the  time  evolution  of  the  critical  surfaces  of  liquid  oxygen  (LOX) 
droplets  under  various  flow  conditions  at  p=100  atm.  Owing  to  the  difference  between 
mass  and  thermal  diffusivity,  the  surface  which  attains  the  critical  mixing  temperature 
usually  regresses  faster  than  that  with  the  critical  mixing  composition.  The  dynamic 
deformation  of  the  droplet  is  substantially  enhanced  by  increasing  the  momentum  carried 
by  the  ambient  flow.  The  effect  of  pressure  on  the  evolution  of  the  critical  surfaces  is 
illustrated  in  Fig.  2.  In  addition  to  enhanced  droplet  deformation,  entrainment  of  gasified 
oxygen  into  the  recirculating  eddies  in  the  wake  region  is  augmented  with  increasing 
pressure.  At  p=400  atm,  the  isocomposition  skirt  even  exhibits  oscillatory  motion.  The 
gaseous  oxygen  entrapped  by  the  recirculating  flow  tends  to  move  forward  and  drives  the 
skirt  to  expand  in  the  cross  steam  direction.  The  convective  flow,  on  the  other  hand, 
suppresses  this  lateral  expansion  and  forces  the  skirt  to  bend  the  stretch  downstream.  As 
a  result,  more  oxygen  is  trapped  into  the  recirculating  eddies,  leading  to  another 
locomotion.  Detailed  discussions  of  this  phenomenon  are  given  in  Refs.  5  and  6. 

2.  Dynamic  Responses  of  Droplet  Vaporization  and  Combustion  to  Ambient  Flow 
Oscillations 

The  dynamic  responses  of  supercritical  droplet  vaporization  and  combustion  to 
ambient  flow  oscillations  have  been  examined  [7-9].  The  analysis  extends  the  droplet 
combustion  model  established  for  steady  environments,  and  accommodates  periodic 
pressure  and  velocity  fluctuations  in  the  ambiance.  The  oscillatory  characteristics  of 
droplet  gasification  and  burning  mechanisms  are  studied  in  detail  over  a  broad  range  of 
pressure.  Results  indicate  that  the  droplet  response  functions  (both  pressure-  and 
velocity-coupled  responses)  increase  progressively  with  pressure  due  to  reduced  thermal 
inertial  at  high  pressures.  Furthermore,  a  rapid  amplification  of  droplet 
vaporization/combustion  response  occurs  when  the  droplet  surface  reaches  its  critical 
mixing  point.  This  phenomenon  may  be  attributed  to  the  strong  sensitivities  of  latent 
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figure  2.  Effect  of  Ambient  Pressure  on  Evolution  of  Droplet  Critical  Sui&ce, 
LOX  /  H2  System  at  J7oo  =  2.5  m/s. 


heat  of  vaporization  and  fluid  transport  properties  to  small  perturbations  near  the  critical 
point. 


In  addition  to  basic  investigation  into  the  thermophysical  processes  involved, 
correlations  of  droplet  response  functions  are  established  in  terms  of  droplet  and 
oscillatory  flow  properties.  This  information  can  be  effectively  used  in  liquid  rocket 
engine  combustion  instability  analysis. 
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The  gasification  and  dynamics  of  an  isolated  liquid  oxygen  (LOX)  droplet  in  a  su¬ 
percritical  hydrogen  stream  has  been  numerically  studied  based  on  the  complete  con¬ 
servation  equations.  The  approach  employs  a  state-of-the-art  treatment  of  thermophys¬ 
ical  properties,  and  accommodates  thermodynamic  nonideality  and  transport  anomaly 
at  high  pressure.  The  analysis  allows  for  a  rigorous  investigation  into  droplet  behav¬ 
ior  over  broad  ranges  of  pressure  and  Reynolds  number.  Detailed  flow  structures  and 
transport  phenomena  are  examined  to  reveal  various  key  mechanisms  underlying  droplet 
gasification  in  a  supercritical,  forced-convective  environment.  In  addition,  correlations 
of  droplet  lifetime  and  drag  coefficient  are  established  in  terms  of  fluid  thermodynamic 
state,  Reynolds  number,  and  gasification  transfer  number. 

Nomenclature 

a  Acceleration  of  droplet  center  of  gravity 

B  Spalding  transfer  nuinber 

Bo  Transfer  number  for  drag  correlation 

Co  Drag  coefficient 

Cp  Constant-pressure  specific  heat 

D  Molecular  diameter 

F  Convective-flux  vectors  in  axial  and  radial  directions,  respectively 
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Ex;,  Fy  DifFusion-fiux  vectors  in  axial  and  radial  directions,  respectively 

e  Specific  internal  energy 

Cf  Specific  total  internal  energy 

T  Scaling  factor  for  transport  property 

/,  h  Scaling  factors  for  corresponding  state  principle 

h  Specific  enthalpy 

kij,  £ij  Binary  interaction  constants  for  energy  and  volume,  respectively 

M  Droplet  mass  or  molecular  weight 

N  Number  of  species 

p  Pressure 

qe  Energy  diffusion  flux 

9m, t  Mass  diffusion  flux  of  species  i 

Q  Dependent  variable  vector 

r  Radial  coordinate 

R  Universal  gas  constant 

Rc  Radius  of  mass-equivalent  sphere 

Re  Reynolds  number 

S  Source  vector 

i  Time 

T  Temperature 

u,  V  Velocities  in  axial  and  radial  directions,  respectively 

Udrop  Velocity  of  droplet  center  of  gravity 

V  Specific  volume 

X  Axial  coordinate 


V 


Molar  volume 
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X  Correction  factor  for  noncorrespondence 

Xi  Mole  fraction  of  species  i 

Yi  Mass  fraction  of  species  i 

Z  Compressibility 

Z  Pseudo-time  variable  vector 

Greek  Symbols 

0  Scaling  factor  for  pseudo-time  variable 

r  Preconditioning  matrix 

AA  Dense-fluid  correction  for  thermal  conductivity 

A;i  Dense-fluid  correction  for  viscosity 

T)  Generalized  coordinates 

9  Shape  factor  for  energy 

A  Thermal  conductivity 

/I  Viscosity 

p  Density 

r  ^  Shear  stress  tensor  or  pseudo  time 

Tf  Droplet  lifetime 

<t>  Shape  factor  for  molecular  size 

u;  Pitzer’s  acentric  factor 

Subscripts 

0  Reference  fluid  or  initial  condition 

big  Large  molecule 


c 


Critical  condition 
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drop 
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small 

00 

Superscripts 

crit 

exc 

T 

* 


Droplet 
Species  i 

Averaged  quantity  of  liquid  droplet 
Mixture 

Reduced  thermodynamic  property 
Droplet  surface 
Small  molecular 
Ambient  flow 


Critical  enhancement 

Dense  fluid  correction 

IVanspose  of  vector 

Reduced  by  value  at  critical  point 

Ideal  gas  property  or  hard  sphere 

Vector  in  generalized  coordinates 
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1.  Introduction 

Supercritical  droplet  gasification  and  combustion  have  long  been  of  serious  concern  in 
the  development  of  high-pressure  combustion  devices  such  as  diesel  engines,  liquid  rocket 
engines,  and  gas  turbines.  Most  of  previous  studies  (e.g.,  Hsieh  ei  al  1991;  Shuen  ei  al 
1992;  Yang  ti  aL  1994;  Lafon  1995;  Lafon  ti  ai  1996)  focused  on  droplet  behavior  at 
quiescent  conditions.  Effects  of  free-stream  velocity  on  droplet  gasification  and  dynamics 
in  forced-convective  environments  have  not  yet  been  addressed  in  detail.  The  present 
work  attempts  to  conduct  a  comprehensive  investigation  into  supercritical  liquid  oxygen 
(LOX)  droplet  gasification  in  a  hydrogen  stream.  The  approach  is  based  on  a  state-of- 
the-art  treatment  of  thermodynamic  and  transport  properties,  and  allows  for  a  rigorous 
analysis  from  first  principles.  Various  fundamental  issues  associated  with  droplet  dy¬ 
namics  and  related  transport  phenomena  at  high  pressures  arc  examined  systematically, 
with  emphasis  placed  on  the  distinct  characteristics  of  cryogenic  fluids. 

In  the  first  part  of  this  study,  a  series  of  calculations  are  carried  out  to  investigate 
LOX  droplet  gasification  in  hydrogen  over  broad  ranges  of  pressure  (100-400  atm)  and 
Re3rnolds  number  (2.5-300).  Detailed  flow  structures  and  transport  phenomena  are  ex¬ 
amined  ta  reveal  various  underl3ring  mechanisms  involved  in  supercritical  droplet  gasi¬ 
fication.  Secondly,  correlatiozis  of  mass  and  momentum  transfer  rates  between  droplet 
and  free  stream  are  established  in  terms  of  fluid  thermodynamic  state,  Reynolds  number, 
and  gasification  transfer  number.  Results  can  be  effectively  used  as  a  physical  submodel 
for  high-pressure  spray  combustion  analyses. 

2.  Theoretical  Formulation 

The  physical  model  of  interest  includes  an  isolated  LOX  droplet  in  a  hydrogen  stream, 
as  shown  in  Fig.  1.  The  initial  droplet  temperature  is  uniformly  distributed  at  a  suberiti- 


6  G.  C.  Hsiao,  V.  Yang  and  J.  S.  Shuen 

cal  value,  but  the  ambient  hydrogen  pressure  and  temperature  are  in  the  thermodynamic 
supercritical  regime  of  oxygen.  Owing  to  the  cryogenic  fluid  property  of  oxygen,  the 
droplet  surface  reaches  its  thermodynamic  critical-mixing  state  almost  instantaneously 
upon  introduction  into  the  hydrogen  flow.  Once  this  occurs,  the  enthalpy  of  vaporiza¬ 
tion  and  surface  tension  vanish.  The  entire  flowfield  (including  both  the  droplet  interior 
and  surrounding  gases)  becomes  essentially  a  continuous  medium,  and  no  well-defined 
liquid/gas  interfacial  boundary  existSf^as^for  a  subcritical  droplet.  The  temperature  and 
density,  as  well  as  their  spatial  gradients,  vary  smoothly  throughout  the  domain  from  the 
liquid  core  to  the  far  field.  As  a  consequence,  a  single-phase  fluid  model  is  used  herein 
to  facilitate  the  analysis. 

The  flowfield  is  assumed  to  be  laminar  and  arisymmetric.  If  body  forces  and  thermal 
radiation  are  ignored,  the  conservation  laws  of  mass,  momentum,  energy,  and  spedes 
concentration  can  be  written  in  the  following  vector  form. 

dQ  ,  d{E-E,)  diF-F,)_  (2.1) 

dt  ^  dx  dr 

The  conservative-variable  vector  Q  is  defined  as 


Q  =  r{p,  pu,  pv,  pet,pYi)'^  (2*2) 

where  superscript  T  denotes  the  transpose  of  vector.  The  physical  properties  p,  (ti,  u), 
and  Yi  are  the  density,  velocity  components  in  the  axial  and  radial  directions,  and  mass 
fraction  of  species  i,  respectively.  The  specific  total  internal  energy  is  written  as 


LOX  DROPLET  GASIFICATION  IN  SUPERCRITICAL  HYDROGEN  FLOWS  7 
where  h  denotes  the  specific  enthalpy  and  p  the  pressure.  The  convective-flux  vectors,  E 

and  F,  in  the  axial  and  radial  directions,  respectively,  take  the  forms 


E  =  r  ^pu,  pu^  +  p,  puw,  (pe,  +  p)u,  puY^ 

(2.4) 

F  =  r  ^pv,  puv,  pv^  +  p,  {pet  +  p)v,  pvY^ 

(2.5) 

The  corresponding  diffusion-flux  vectors  are 


F.  =  r(0  ^f’xx  "I"  ^^xr  "t"  (9«)*i  (9m, «)®^  (^•®) 

Fv  —  f(0,TxrtTrr,UTgf  +  VTff  +  (9e)rt  (9m,«)r^  (2*'^) 

where  the  normal  stresses  {txi  and  r^r)  and  shear  stresses  (r^r  and  Vn)  are 


dux  2 

=  2^:^ -3^ 


dux  ,1  d  .  . 

IT 


dur  2  r  dux  ,  1  ^  /  0 

r< 

=  Trx  =  fi\- 


\dUx 

dr 


dUr 

dx 


(2.8) 

(2.9) 

(2.10) 


The  energy  and  mass  diffusion  terms  (g.  and  9m, t)  calculated  using  Fourier’s  and 

Pick’s  laws,  respectively.  The  source  vector  S  results  from  the  axisymmetric  geometry 
and  can  be  written  as 


S  = 


,0,0)' 


(2.11) 


where 


rt  2 

T90  =2p— 

r  0 


dUr  ,15,  . 

aT 


(2.12) 
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3.  Property  Evaluation 

One  major  challenge  in  the  analysis  of  supercritical  droplet  behavior  is  the  establish¬ 
ment  of  a  unified  property  evaluation  scheme  capable  of  treating  thermophysical  prop¬ 
erties  of  pure  substances  and  muctures  over  the  entire  fluid  thermodynamic  state  &om 
compressed  liquid  to  dilute  gas.  At  high  pressures,  models  normally  used  to  represent 
ideal-gas  behavior  may  encounter  significant  dfficulties.  From  the  microscopic  point  of 
view,  the  intcrmolecular  mean  free  paths  tend  to  decrease  with  increasing  pressure,  and 
consequently  the  molecular  volume  and  intermolecular  forces  are  no  longer  negligible  as 
in  the  case  for  idealized  fluids.  For  convenience,  each  property  may  be  expressed  as  the 
sum  of  an  ideal-gas  property  at  the  same  temperature  and  a  thermodynamic  departure 
function  which  takes  into  account  the  dense-fluid  correction. 

In  the  current  study,  a  modified  Benedict- Webb- Rubin  (BWR)  equation  of  state  pro¬ 
posed  by  Jacobsen  and  Stewart  (1973)  is  used  to  represent  the  fluid  p-V-T  behavior  due 
to  its  superior  performance  over  conventional  cubic  equations  of  state  (Prausnitz  ei  al. 
1986).  This  equation  of  state  has  been  extremely  valuable  in  correlating  both  liquid  and 
vapor  thermodynamic  and  volumetric  data;  however,  the  temperature  constants  involved 
are  available  only  for  a  limited  number  of  pure  substances  (Cooper  ei  ai  1967;  Orye 
1969;  Perry  ei  al,  1984).  To  overcome  this  constraint,  an  extended  corresponding-state 
(ECS)  principle  (Ely  &  Hanley  1981,  1983)  is  used  herein.  The  basic  idea  is  to  assume 
that  the  properties  of  a  single-phase  fluid  can  be  evaluated  via  conformal  mappings  of 
temperature  and  density  to  those  of  a  given  reference  fluid.  As  a  result,  only  the  BWR 
constants  for  the  reference  fluid  are  needed.  For  a  multicomponent  system,  accounting 
for  changes  in  properties  due  to  mixing  is  much  more  complicated.  A  pseudo  pure- 
substance  model  is  adopted  to  evaluate  the  properties  of  a  mixture,  treating  the  mixture 
as  a  single-phase  pure  substance  with  its  own  set  of  properties  evaluated  via  the  ECS 
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principle.  This  method  improves  prediction  accuracy  and  requires  only  limited  data  (i.e., 
critical  properties  and  Pitzer’s  acentric  factor)  for  each  constituent  component. 

Successful  application  of  the  corresponding-state  argument  for  the  evaluation  of  fluid  p- 
V-T  properties  also  encourages  similar  improvement  in  the  prediction  of  thermophysical 
data.  In  the  following,  a  brief  summary  of  the  corresponding-state  method  in  conjunction 
with  the  mixture  combining  rule  is  first  given,  followed  by  the  B  WR  equation  of  state  for 
the  reference  fluid.  Finally,  the  techniques  for  evaluating  thermodynamic  and  transport 
properties  axe  addressed. 


3.1.  Extended  Corresponding-Siaie  Principle 

The  extended  corresponding-state  model  of  Ely  and  Hanley  (1981, 1983)  is  used  to  evalu¬ 
ate  volumetric  and  transport  properties  of  a  mixture  over  its  entire  thermodynamic  fluid 
state.  The  scheme  assumes  that  the  configurational  properties  (such  as  temperature, 
density,  viscosity,  thermal  conductivity,  etc.)  of  a  single-phase  mixture  can  be  equated 
to  those  of  a  hypothetical  pure  fluid,  which  are  then  evaluated  via  corresponding  states 
with  respect  to  a  given  reference  fluid.  For  example,  the  viscosity  of  a  mixture,  fimy  can 
be  related  to  that  of  a  reference  fluid,  /iq,  at  the  corresponding  thermodynamic  state  as 

Mm  (p,  T*)  =  Mo(po  yTo)  Pfi  (3.1) 

where  represents  the  mapping  function.  The  correspondence  of  temperature  and 
density  between  the  mixture  of  interest  and  the  reference  fluid  (denoted  respectively  by 
subscripts  m  and  O)  can  be  characterized  by  the  following  two  scaling  factors. 

=  ^  ;  fim  =  —  (3.2) 

-to  p 

The  former  represents  the  conformation  of  potential  distribution  of  energy,  while  the 
latter  characterizes  the  effect  of  mixture  molecular  size.  Assuming  that  all  components 
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in  a  mixture  obey  the  ECS  principle,  the  mixing  rules  for  a  multicomponent  system  can 
be  expressed  as 

N  N 

»  J 

N  N 

XiXjfijhij  (3.3) 

»  J 

where  N  is  the  number  of  species  and  Xi  the  mole  fraction  of  species  i.  The  two  binary 
scaling  factors,  fij  and  hij,  are  defined  as 

fij  =  (/.7;)‘^'(1  -  kij) 

=  +  (3.4) 

where  kij  and  fy  are  binary  interaction  parameters  reflecting  the  effects  of  energy  and 
molecular  size,  respectively.  In  order  to  model  the  quantum  behavior  of  hydrogen,  inter¬ 
action  parameters  for  a  hydrogen-contMning  mixture  are  specially  treated  by  a  general¬ 
ized  correlation  proposed  by  Valderrama  and  Reyes  (1983). 

kij  =  A-  B/TJ 

iij  =  l  +  0.071n[(l>4,j/D,„aH)]  (3,5) 

and 

[Dus/D,„,u]  = 


A  =  0.1805  +  Z.21uj  +  2.437W? 

B  =  0.1323  +  0.5507wj-  +  3.5994w? 


(3.6) 
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where  \DuglD,mau]  is  always  greater  than  1  and  Vc  is  the  critical  molar  volxime.  The 

variable  T*  is  the  reduced  temperature  and  u  Pitzer’s  acentric  factor.  The  subscripts  t 
and  j  denote  the  hydrogen  and  non-hydrogen  components,  respectively. 

The  scaling  factors  (/j  and  hi)  for  each  individual  component  can  be  obtained  by  a 
two-parameter  corresponding  state  principle  (Ely  &  Hanley  1981,  1983). 

fi  =  (rc..7r,,o)  ei{T:,viwi) 

fti  =  (Ve.</Vc.o)^.(ir,V:,a;.)  (3.7) 

where  T*  is  the  critical  temperature,  and  Oi  and  are  the  so-called  shape  factors  which 
are  functions  of  Pitzer’s  acentric  factor  (w.  )  and  reduced  temperature  and  molar  volume 

)•  analytical  expressions  for  9i  and  proposed  by  Leach  et  al.  (1968) 
take  the  forms 

«i(ir,v;,u;.)  =  1  +  (Wi  -  wo)  Fijr.v;) 

<f>i{T^ ,  V*, Wj)  =  [1  +  (wj  -  Wo)  G{'I^ , Ni)]Ze,o/Ze,i  (3.8) 

where  Z  is  the  compressibility  factor.  The  variables  F  and  G  are 

V,-)  =  oi  -I-  6i  lni;^(ci  -t-  di/i;^)(V,^  -  0.5) 

G(77 .  V,’)  =  a2(V/'  -I-  62)  +  C2(Vt  -I-  d2)  In  IT*"  (3.9) 

=  mm{2,mox{IJ*,0.5}} 

Vt  =  mm{2,  max{V,r,  0.5}}  (3.IO) 

where  constants  Oj,  6,-,  a,  and  d,-  are  coefficients  for  shape-factor  correlations. 
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3.2.  Equation  of  State 


Under  the  assumption  of  the  ECS  principle,  the  density  of  a  mixture  can  be  evaluated 
by 


where  po,  To,  and  po  denote  the  corresponding  density,  temperature,  and  pressure  of  the 
reference  fluid,  respectively.  Since  the  temperature  at  the  conformal  state  is  by 

Eq.  (3.2),  the  corresponding  pressure  can  be  derived  based  on  the  general  compressibility 
theory  (i.e.,  Z=Z(T*,  p*)), 


P.  =  p(^)  (M2) 

To  ensure  the  accuracy  of  the  density  prediction,  a  generalised  BWR  equation  of  state 
(Jacobsen  &  Stewart  1973)  is  adopted  for  the  reference  fluid. 

9  15 

Po{T,p)  =  ^  a„(T’)p”  +  ^  0„(T)p2'‘-^’’e“'''’*  (3.13) 

n=l  n=10 

where  y  is  0.04,  and  the  temperature  coefficients  ai(T)  depend  on  the  reference  fluid  of 
concern. 

Although  this  equation  of  state  must  be  solved  iteratively  for  density  at  given  pressure 
and  temperature,  the  prediction  covers  wide  ranges  of  thermodynamic  states,  and  as  such 
promotes  the  establishment  of  a  unified  evaluation  scheme  of  thermophysical  properties. 
Figure  2  shows  the  comparison  of  oxygen  density  between  experimental  data  (Sychev  ei 
al  1987)  and  the  prediction  by  the  BWR  equation  of  state  in  conjunction  with  the  EGS 
principle.  The  reference  fluid  is  selected  as  propane  herein  due  to  the  availability  of  suffi¬ 
ciently  reliable  data  correlated  over  a  wide  range  of  experimental  conditions.  The  result 
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shows  excellent  agreement  over  the  entire  fluid  state,  from  compressed  liquid  to  dilute 

gas.  Figure  3  presents  the  relative  errors  of  density  prediction  based  on  three  commonly 
used  equations  of  state,  namely,  the  Benedict- Webb-Rubin  (BWR),  Peng-Robinson  (PR), 
and  Soave-Redlich-Kwong  (SRK).  The  ECS  principle  is  embedded  into  the  evaluation 
procedure  of  the  BWR  equation  of  state  and  shows  its  superior  performance  with  the 
maximum  relative  error  at  1.5%  for  the  pressure  and  temperature  ranges  under  consid¬ 
eration.  On  the  other  hand,  the  SRK  and  PR  equations  of  state  yield  TnAviTnum  errors 
aroimd  13%  and  17%,  respectively. 


3.3.  Thermodynamic  Properties 

Thermod3mamic  properties  are  of  vital  importance  in  the  analysis  of  supercritical  droplet 
behavior.  Variations  in  these  properties  can  often  be  related  to  local  fluid  properties 
of  pressure,  temperature,  and  species  concentration.  Properties  like  enthalpy,  internal 
energy,  and  specific  heat  of  a  mixture  are  expressible  as  the  sum  of  the  ideal-gas  property 
and  a  correction  term  accounting  for  thermodynamic  nonideality. 


(3.14) 


(3.15) 


(3.16) 


where  superscripts  0  refer  to  ideal-gas  properties.  The  second  terms  on  the  right  sides  of 
Eqs.  (3.14)-(3.16)  denote  the  thermodynamic  departure  functions,  and  can  be  obtained 
from  the  equation  of  state  described  previously. 
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3.4.  Transpori  Pro’periits 


Estimation  of  viscosity  and  thermal  conductivity  can  be  made  by  means  of  the 
corresponding-state  method  described  in  Sec.  3.1.  The  corresponding-state  argument 
for  the  viscosity  of  a  mixture  can  be  written  in  its  most  general  form  as 

=/1o(POj2o)  Tfx  X/4  (^•17) 

The  quantity  is  a  correction  factor  accounting  for  the  effect  of  noncorrespondence. 
Its  magnitude  is  always  close  to  unity,  based  on  the  modified  Enskog  theory  (Ely  1981). 
The  scaling  factor  takes  the  form 


—  (Affn/<^o)^  /m  ^  (3.18) 

where  Mo  is  the  molecular  weight  of  the  reference  fluid,  and  the  value  of  the  mixture 
denoted  by  Mm  can  be  evaluated  by 

Mm  =  (3.19) 

*-1=1  i  =  l  J 

The  scaling  factors  fm^  /ij»  and  fijj  follow  the  same  definitions  in  Eqs.  (3.2)  and 
(3.4).  A  binary  mixing  rule  is  used  for  the  molecular  weight. 


Mi^  =  2M,Mj/(M<  -h  Mj)  (3.20) 

where  M,*  is  the  molecular  weight  of  component  i. 

Because  of  the  lack  of  a  complete  molecular  theory  for  describing  transport  properties 
over  a  broad  regime  of  fluid  phases,  it  is  generally  accepted  that  viscosity  and  thermal 
conductivity  be  divided  into  three  contributions  and  correlated  in  terms  of  density  and 
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temperature  (Vesovic  &  Wakeham  1991).  For  instance,  the  viscosity  of  the  reference  fluid 

is  written  as  follows. 

Po{po,To)  =  Po{To)  +  To)  +  Afiy‘**{po,To)  (3.21) 

The  first  term  on  the  right-hand  side  represents  the  value  in  the  dilute-gas  limit,  which  is 
independent  of  density  and  can  be  accurately  predicted  by  kinetic-theory  equations  (Ely 
&  Hanley  1981,  1983).  The  second  term  is  the  excess  viscosity  which,  with  the  exclusion 
of  unusual  variations  near  the  critical  point,  characterizes  the  deviation  from  po  for  a 
dense  fluid.  The  third  term  refers  to  the  critical  enhancement  which  accounts  for  the 
anomalous  increase  above  the  background  viscosity  (i.e.,  the  sum  of  /ig  and  A/i§*®)  as 
the  critical  point  is  approached.  However,  the  theory  of  nonclassical  critical  behavior 
predicts  that,  in  general,  properties  that  diverge  strongly  in  pure  fluids  near  the  critical 
points  will  diverge  only  weakly  in  mixtures  due  to  the  different  physical  criteria  for 
criticality  in  a  pure  fluid  and  a  mixture  (Levelt  Sengers  1991).  Because  the  effect  of 
critical  enhancement  is  not  well-defined  for  a  mixture  and  is  likely  to  be  small,  the  third 
term  is  not  considered  in  the  existing  analyses  of  supercritical  droplet  gasification. 

Evaluation  of  thermal  conductivity  must  be  carefully  conducted  for  two  reasons:  (1) 
the  one-fluid  model  must  ignore  the  contribution  of  diffusion  to  conductivity,  and  (2)  the 
effect  of  internal  degrees  of  ffeedom  on  thermal  conductivity  cannot  be  correctly  taken 
into  account  by  the  corresponding-state  argument.  As  a  result,  thermal  conductivity  of 
a  pure  substance  or  mixture  is  generally  divided  into  two  contributions  (Ely  &  Hanley 
1983), 


A„.(p,r)  =  A'„(r)-hA"(p,r) 


(3.22) 


The  former,  arises  from  transfer  of  energy  via  the  internal  degrees  of  freedom, 
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while  the  latter,  A"  (p,r),  is  due  to  the  effects  of  molecular  collision  or  translation.  For 


a  mixture,  A(^(T)  can  be  evaluated  by  a  semi-empirical  mixing  rule. 


=  E  E  (3.23) 

*  J 

where  AJ^  is  a  binary  thermal  conductivity  defined  as 

=  2[a:.(T)-'  +  AJ(T)-1]  (3.24) 

The  internal  contribution  of  component  i  is  calculated  by  the  modified  Eucken  correlation 
(Ely  k  Hanley  1983)  for  polyatomic  gases  as 


A;.(r)  =  i.32^(c°i-|ie) 


(3.25) 


where  the  properties,  i?,  and  Cp  i,  are  the  dilute-gas  viscosity,  gas  constant,  and 
ideal-gas  heat  capacity  of  component  i,  respectively.  The  collisional/translational  con¬ 
tribution,  A"  (p,  r),  is  evaluated  via  the  extended  corresponding-state  method. 


^m{Pi T)  =  Ao(po,  To)  Tx  XX  (3.26) 

where  is  the  scaling  factor  identical  to  that  used  in  Eq.  (3.18)  for  mixture  viscosity. 
A  small  correction  factor  xx  is  used  to  compensate  the  noncorrespondent  behavior  of 
the  mixture,  and  can  be  estimated  using  the  modified  Enskog  theory  and  fluid  p-V-T 
properties.  The  quantity  Ao(po,To)  is  the  equivalent  contribution  for  the  reference  fluid 
at  the  corresponding  state,  which  can  be  further  divided  into  three  parts: 


KiPo.  To)  =  A;(r)  -h  AAg-(po ,  To)  +  AA-**(po,  To)  (3.27) 


where  A;(r)  is  the  value  in  the  dilute-gas  limit,  AAg®^  the  excess  thermal  conductivity, 
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«d  AAS-  .he  .riMc,  ^ 

.l.»..l  conductivU,  .pp^  „  b. 

for  critical  .nl„..c«n«.t,  AAS«',  b  act  iachdcd  ia  the  p,«„ai 
E-ttaarioa  ,f  tka  ^  ^  ^  ^  ^  ^ 

«».  chaUca^g  taak  .kaa  .valaatias  tk.  ,tk„  riaaepari  p„p„i„,  <,a.  t.  tk.  lack  ofa 
fa^al  tk«,„  ac  „.a  a  «.ca„.lcall,  baaed  cattriatlaa.  makaaki  (1974)  haa  aaaaried 
a  aiaapl.  ackaat.  fa,  pradicria,  tha  biaaq,  diHarivit,  af  a  daaa.  Bald  b,  maaaa 
ttf  a  ca„«,p,adiag^U.a  app,«,ch.  Thia  ach«.,  appaar.  ,a  b.  th.  mat  c„pi„a  to 
d«.,  and  ha.  demaaatoatod  ntodatoto  .«„,„  in  toe  lhaitod  aantoe,  af  teak,  caadnetod. 
The  .ppa>ach  ptoc«»ia  ia  toa  atop,.  Pleat,  th.  biato,  atoa,  diilaririt,  af  a  dilute  ,a. 

1.  aktala^i  todag  toe  Cl.apatoa-Ea.kag  tkeat,  ia  caapa.c«aa  trith  th.  iata,-a.al«ari.. 

i»to..Ul  fkactlaa.  Tk.  calculatad  data  i.  toaa  ca^ctad  ia  actoHaac.  trith  a  gaaatoliaki 
chart  in  terms  of  reduced  temperature  and  pressure. 

4.  Numerical  Algorithm 

Depict  toparitotiaa  aad  ca.abu.tiaa  tjpically  iavala.  f  aid  atatiaua  ia  a  toladty  toag. 
lacula,  dilfua.aa  to  lew  aubaaaic  speed.  Caatontparary  nuntorical  algarithma 
devriapril  fa,  cutoptoadbl.  Baw.  a,,  aftoa  iaaffactiva  a,  auri.  a  la.-yalacity  caadiriaa. 
Thcto  „a  ,dl.,w»,priawl  toaaaaa  fa,  thia  difflculty  (Me,kl.  b  Obai  1985).  Pirn,  th. 
rigeawriuc.  af  th.  ayataa.  bacaaa.  ati*  a.  la.  Ba,  veladti.,,  toetaby  advatoaly  aifactiag 
th.  caatotgaac.  ckaeactariatic.  af  auataricri  «>l„Haua.  Sacaad.  to.  ptoaauto  k. 

th.  a.am.atoa.  wpariiaa.  hacaa..  riagula,  a.  th.  Mtoh  auatba,  H.p,«toha.  and 

consequently  result  in  a  large  round-off  error. 

Ia  aalat  to  cacumvaat  to.  ahava  camputatiaaal  difficultiaa,  a  faUy  iaiplicit  dual  tin,.- 
riappbg  ia..g,«laa  atatoad  (Withiagtaa  a,  C.  1991.  Shaaa  „  .1  1992)  i.  «,.ptod 
ia  tha  p,««.t  study.  Th.  acluaaa  1,  aatobllahed  ia  tw,  stop.  Pi,.., 
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terms  are  used  in  the  momentum  equations  to  minimize  round-off  errors  at  low  Mach 
numbers.  The  second  step  incorporates  a  set  of  well-conditioned  artificial  terms  into 
the  conservation  laws  to  enhance  numerical  efficiency  and  stability.  An  implicit  iterative 
procedure  follows  to  achieve  a  converged  solution  in  the  pseudo-time  domain,  which 
corresponds  to  a  temporal-accurate  solution  in  the  physical-time  domain.  The  resulting 
conservation  laws  in  generalized  coordinates  ,  rj)  are  written  as 


,  d{F-F,)  . 

dr  ^  dt^  ^  -  * 

where  r  represents  the  pseudo  time.  The  vectors  Q,  E,  E„,  F,  and  F„ 


(4.1) 

are  defined  as 


E=j{^tQ  +  ^sE  +  ^rF) 
Ev  =  •j{^xEv 
F  =  +  VxE  -I-  ffrF) 

Fv  —  -j{r)sEv  -t-  VrFv) 


(4.2) 


where  J  is  the  Jacobian  matrix  for  coordinate  transformation.  Variables  and  »?<  are 
the  grid-speed  terms,  and  jj®,  and  the  metric  terms.  The  primitive  variable 

vector,  Z,  the  preconditioning  matrix,  T,  and  the  source  term,  5,  are  expressed  as 
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where  pg  and  0  are  the  gauge  pressure  and  its  scaling  factor,  respectively. 

One  significant  advantage  of  the  dual  time-stepping  integration  method  is  that  the 
convergence  of  the  iterative  process  is  determined  by  the  eigenvalue  characteristics  in  the 
pseudo-time  domain,  not  by  the  original  eigenvalues  which  become  disparate  at  low  flow 
velocity.  This  feature  allows  flenbility  in  the  selection  of  time  step  sizes  in  both  time 
frames.  The  physical-time  step  is  determined  based  on  the  characteristic  evolution  of  the 
unsteady  flow  under  consideration,  while  the  pseud<vtime  step  depends  on  the  numerical 
stability  of  the  algorithm  and  can  be  adjusted  to  obtain  optimum  convergence. 

5.  Results  and  Discussions 

The  model  described  in  the  preceding  sections  has  been  applied  to  study  the  gAnifiration 
characteristics  of  liquid  oxygen  (LOX)  droplets  in  supercritical  hydrogen  streams.  The 
ambient  hydrogen  temperature  is  taken  to  be  1000  K.  The  initial  droplet  diameter  is 
100  pm,  and  the  initial  droplet  temperature  is  100  K,  which  corresponds  to  the  LOX 
injection  temperature  in  many  operational  rocket  engines  such  as  the  Space  Shuttle 
main  engines  and  the  Vulcain  engines  used  in  the  Ariane  5  launch  vehicles.  In  order 
to  examine  the  effects  of  ambient  flow  conditions  on  droplet  behavior  and  aag/winf^ 
interphase  transport,  a  parametric  study  covering  a  wide  range  of  pressure  (100-400  atm) 
and  velocity  (0.1-15  m/s)  is  conducted.  The  solid  dots  in  Fig.  4  mark  the  scenarios  under 
consideration,  with  the  corresponding  free-stream  Reynolds  number  based  on  the  initial 
droplet  diameter  {Re  =  U„Do/voo)  shown  in  the  y-axis. 

Within  the  pressure  range  of  the  current  interest,  the  droplet  initial  temperature  is 
slightly  lower  than  the  critical  mixing  temperature  of  the  oxygen/hydrogen  system  (i.e., 
Te  =  142.79  /f  at  p  =  100  atm).  Thermodjmamic  criticality  is  reached  at  the  droplet 
surface  almost  instantaneously  upon  its  introduction  to  the  flow  due  to  rapid  heat  transfer 
from  the  ambient  flow  (Yang  et  al.  1994,  Lafon  et  ai  1996).  Once  this  occurs,  both  the 
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liquid  and  gas  phases  coexist  on  the  droplet  surface  at  the  critical  mixing  condition.  The 
enthalpy  of  vaporization  and  surface  tension  reduce  to  zero,  and  the  mixture  properties 
as  well  as  their  spatial  gradients  vary  smoothly  across  the  surface.  However,  the  droplet 
interior  still  remains  at  the  liquid  state  with  a  subcritical  temperature  distribution.  The 
subsequent  gasification  process  becomes  difFusion/convection  controlled,  and  the  effect 
of  interfacial  thermodynamics  diminishes. 

Several  snapshots  of  the  flowfields  under  various  conditions  are  presented  in  Figs.  5- 
7.  Figure  5  shows  six  frames  of  isotherms  and  isopleths  of  oxygen  concentration  at  a 
convective  velocity  of  2.5  m/s  and  an  ambient  pressure  of  100  atm.  The  free  stream 
Reynolds  number  Re  is  30.56  based  on  the  initial  droplet  diameter.  Soon  after  the 
introduction  of  the  droplet  into  the  hydrogen  stream,  the  flow  adjusts  to  form  a  boundary 
layer  near  the  surface.  The  gasified  oxygen  can  not  accumulate  on  the  front  surface  of  the 
droplet,  and  is  carried  downstream  through  convection  and  mass  diffusion.  The  evolution 
of  the  temperature  field  exhibits  features  distinct  from  that  of  the  concentration  field  due 
to  the  disparate  time  scales  associated  with  thermal  and  mass  diffusion  processes  (i<e., 
Lewis  number  ^1).  The  thermal  wave  penetrates  into  the  droplet  interior  with  a  pace 
faster  than  does  the  surrounding  hydrogen  species.  Since  the  liquid  core  possesses  large 
momentum  inertia  and  moves  slower  compared  with  the  gasified  oxygen,  at  t  =  0.79  ms, 
the  droplet  (delineated  by  the  dark  region  in  the  temperature  contours)  reveals  an  olive 
shape  while  the  oxygen  concentration  contours  deform  into  a  crescent  shape  with  the  edge 
bent  to  the  streamwise  direction.  At  f  =  1.08  ms,  the  subcritical  liquid  core  disappears, 
leaving  behind  a  puff  of  dense  oxygen  fluid  which  is  convected  further  downstream  with 
increasing  velocity  until  it  reaches  the  momentum  equilibrium  with  the  ambient  hydrogen 
flow. 


Figure  6  presents  a  sequence  of  droplet  evolution  as  does  Fig.  5,  but  at  an  ambient 
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pressure  of  400  atm.  The  effect  of  pressure  on  droplet  dynamics  can  be  determined  by  a 

direct  comparison  between  these  two  figures.  Several  distinct  phenomena  are  observed. 
First,  the  droplet  deformation  is  enhanced  by  increasing  the  ambient  pressure.  Here, 
the  corresponding  Reynolds  number  is  3.75  times  larger  than  that  at  100  atm.  The 
hydrogen  stream  carries  higher  momentum  and  consequently  exerts  stronger  force  on  the 
droplet,  and  as  such  promotes  its  deformation.  Second,  the  droplet  lifetime  decreases 
with  increasing  ambient  pressure  because  the  higher  convective  momentum  increases  the 
contact  surface  exposed  to  the  hot  stream,  thereby  facilitating  the  transfer  of  thermal 
energy.  Third,  the  increased  ambient  pressure  advances  the  gradients  of  temperature  «-nd 
oxygen  concentration.  Although  the  blowing  effect  due  to  droplet  gasification  is  stronger 
in  a  high-pressure  environment,  the  large  momentum  transfer  reduces  the  convective 
resident  time  and  forms  a  thinner  boundary  layer  near  the  surface. 

Figure  7  depicts  a  scenario  identical  to  that  of  Fig.  5,  but  with  an  increased  convective 
velocity  of  15  m/s.  The  droplet  evolution  is  substantially  different,  exhibiting  several 
distinct  modes  such  as  deformation,  viscous  stripping,  and  breakup.  When  the  droplet  is 
introduced  into  this  strong  convective  flow  {Re  =  200),  fast  vorticity  generation  caused  by 
the  velocity  difference  between  the  liquid  core  and  ambient  flow  promotes  the  formation 
of  an  attached  eddy  behind  the  droplet  and  induces  a  slightly  higher  pressure  region. 
Meanwhile,  the  free  stream  exerts  dynamic  loading  continuously  to  the  front  surface. 
As  a  result,  the  droplet  deforms  and  extends  in  the  direction  normal  to  the  external 
flow.  The  flattened  edge  then  turns  to  the  streamwise  direction,  stretches  downstream, 
and  forms  a  skirt.  The  length  of  the  skirt  increases  with  time,  and  is  dictated  by  the 
momentum  transfer  with  the  ambient  flow.  At  t=92  /is,  an  of  gaseous  oxygen 

is  detached  from  the  growing  skirt  due  to  the  local  flow  motion  and  volume  dilatation 
resulting  from  the  penetration  of  thermal  wave  into  the  skirt.  The  stripped  oxygen. 
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entrained  by  the  recirculating  flow,  further  bends  toward  the  rear  surface  of  the  droplet. 

It  carries  the  momentum  from  the  attached  eddy  and  consequently  causes  the  droplet 
to  deform  into  a  spherical-cape  or  crescent  shape.  The  associated  cusp-like  rim  may 
catalyze  the  vorticity  generation  and  further  enhance  the  strength  of  the  recirculating 
eddy.  Since  the  viscosity  of  the  droplet  interior  is  greatly  reduced  by  the  energy  transfer 
from  the  hot  stream,  at  certain  stage  during  the  droplet  lifetime,  the  liquid  oxygen  breaks 
up  and  forms  a  secondary  ring  around  its  core.  Another  significant  phenomenon  is  the 
distribution  of  oxygen  concentration.  The  gasified  oxygen  stripped  from  the  edge  of  the 
droplet  leaves  the  recirculating  region  and  is  swept  downstream,  forming  a  cylindrical 
tail. 

Figures  8(a)-(d)  summarize  the  streamline  patterns  and  oxygen  concentration  contours 
of  the  four  different  modes  commonly  observed  in  supercritical  droplet  gasification.  The 
droplet  may  either  remains  in  a  spherical  configuration,  deforms  to  an  olive  shape,  or  even 
breaks  up  into  fragments,  depending  on  the  local  flow  conditions.  Unlike  low-pressure 
cases  in  which  the  large  shear  stress  at  the  gas-liquid  interface  induces  internal  flow  cir¬ 
culation  in  the  liquid  core  (Prakash  Sc  Sirignano  1978,1980),  no  discernible  recirculation 
takes  place  in  the  droplet  interior,  regardless  of  the  Reynolds  number  and  deformation 
mode.  This  may  be  attributed  to  the  diminishment  of  surface  tension  at  supercritical 
conditions.  In  addition,  the  droplet  regresses  so  fast  that  a  fluid  element  in  the  interphase 
region  may  not  acquire  the  time  suflicent  for  establishing  an  internal  vortical  flow  before 
it  gasifies.  The  rapid  deformation  of  the  droplet  configuration  further  precludes  the  ex¬ 
istence  of  stable  shear  stress  in  the  liquid  core  and  consequently  obstructs  the  formation 
of  recirculation. 

The  spherical  mode  shown  in  Fig.  8(a)  typically  occurs  at  very  low  Reynolds  numbers. 
Although  flow  separation  is  encouraged  by  LOX  gasification,  no  recirculating  eddy  is 
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found  in  the  w&ke  behind  the  droplet.  The  vorticity  generated  is  too  weak  to  form  any 

confined  eddy.  When  the  ambient  velocity  increases  to  1.5  m/s,  the  droplet  deforms  into 
an  olive  shape  with  a  recirculating  ring  attached  behind  it,  as  shown  in  Fig.  8(b).  Owing 
to  the  droplet  deformation  and  gasification,  the  threshold  Re3molds  number  above  which 
the  recirculating  eddy  forms  is  considerably  lower  than  that  for  a  hard  sphere.  The  latter 
case  requires  a  minimum  Reynolds  number  of  20  based  on  the  numerical  results  of  LeClair 
et  al.  (1970)  and  Dennis  e<  al.  (1971).  Figure  8(c)  depicts  the  flow  structure  with  viscous 
stripping  at  an  ambient  velocity  of  5  m/s,  showing  an  oblate  droplet  with  a  stretched 
vortex  ting.  The  flattened  edge  of  the  droplet  enhances  the  strength  of  the  recirculating 
*ddie8  and  as  such  increases  the  viscous  shear  stress  dramatically.  Consequently,  a  tViin 
sheet  of  fluid  is  stripped  off  from  the  edge  of  the  droplet  and  swept  toward  the  outer 
boundary  of  the  recirculating  eddy.  At  a  very  high  ambient  velocity  of  15  m/s,  droplet 
breakup  takes  place,  as  clearly  shown  in  Fig.  8(d).  The  hydrogen  flow  penetrates  through 
the  liquid  phase,  and  divides  the  droplet  into  two  parts:  the  core  disk  and  surrounding 
ring.  The  vortical  structure  in  the  wake  region  expands  substantially  as  a  result  of  the 
strong  shear  stress. 

Owing  to  the  difference  between  mass  and  thermal  diffusivity,  the  critical  Tni-ring  state 
cannot  be  sustained  on  the  droplet  surface  after  the  occurrence  of  thermodynamic  criti¬ 
cality.  The  surface  which  attuns  the  critical  mixing  temperature  usually  regresses  faster 
than  that  with  the  critical  mixing  composition.  Fig.  9  shows  the  temporal  evolution 
of  these  two  types  of  critical  surfaces  under  various  flow  conditions  at  p  =  100  atm. 
The  solid  lines  represent  the  instantaneous  isotherms  of  the  critical  mixing  temperature, 
while  the  dashed  lines  mark  the  surfaces  of  the  critical  mixing  composition.  The  val¬ 
ues  of  critical  mixing  properties  can  be  obtained  through  a  phase  equilibrium  analysis 
(Yang  ei  al.  1994)  and  are  given  in  Table  1.  The  dynamic  deformation  of  the  droplet 
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Table  1.  Critical  Mixing  Properties  oi02/H2  System. 


pressure 

(atm) 

critical  mbdng 
temperature  (K) 

critical  mixing 
composition  of  O2 

100 

142.79 

0.735 

200 

127.21 

0.561 

400 

116.0 

0.496 

is  substantially  enhanced  by  increasing  the  momentum  carried  by  the  ambient  flow.  For 
example,  at  one  half  of  the  droplet  lifetime,  the  droplet  aspect  ratio  based  on  the  surface 
of  the  critical  mixing  temperature  (defined  as  the  ratio  of  the  length  in  the  vertical  axis 
to  the  thickness  in  the  axis  of  symmetry)  increases  from  2.1  at  Uoo  =  2.5  m/s  to  3.2  at 
C/qo  =  5  m/ s.  The  effect  of  pressure  on  the  evolution  of  the  critical  surfaces  is  illustrated 
^  addition  to  enhanced  droplet  deformation,  entrainment  of  gasified  oxygen 

into  the  recirculating  eddies  in  the  wake  region  is  augmented  with  increasing  pressure. 
At  p  :=  400  atm,  the  isocomposition  skirt  even  exhibits  oscillatory  motion.  The  gaseous 
oxygen  entrapped  by  the  recirculating  flow  tends  to  move  forward  and  drives  the  skirt  to 
expand  in  the  cross  stream  direction.  The  convective  flow,  on  the  other  hand,  suppresses 
this  lateral  expansion  and  forces  the  skirt  to  bend  and  stretch  downstream.  As  a  result, 
more  oxygen  is  trapped  into  the  recirculating  eddies,  leading  to  another  locomotion. 

Figure  11  presents  the  effect  of  free-stream  velocity  on  the  instantaneous  variation  of 
droplet  residual  mass  at  p=100  atm,  defined  herein  as  the  mass  confined  by  the  isother- 
mal  surface  at  the  critical  mixing  temperature.  In  a  slowly  convective  stream  (e.g., 
Uoo  =0.2  m/s),  heat  conduction  plays  an  important  role  in  determining  the  gasification 
characteristics.  The  rapid  gasification  in  the  early  stage  of  the  droplet  lifetime  results 
from  the  large  temperature  gradient  near  the  liquid/gas  interphase.  The  gasified  oxygen 
then  diffuses  outward  and  hinders  the  penetration  of  thermal  wave  into  the  droplet,  con¬ 
sequently  rendering  a  slower  gasification  rate.  At  high  Reynolds  numbers,  enhanced  heat 
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ajid  mass  transfer,  along  with  severe  droplet  deformation  caused  by  the  strong  dynamic 

loading  from  the  approaching  flow,  considerably  facilitates  the  gasification  process.  The 
droplet  disappears  within  a  much  shorter  period  of  time.  Figure  12  shows  the  effect 
of  pressure  on  droplet  gasification  behavior.  As  the  pressure  increases,  the  convective 
heat  transfer  also  increases  since  the  fieestream  Reynolds  number  (/Ee  =  pU L Jp)  varies 
almost  linearly  with  pressure,  the  viscosity  being  a  weak  function  of  pressure.  In  addi¬ 
tion,  the  critical  mixing  temperature  used  to  define  the  droplet  surface  decreases  with 
increasing  pressure.  These  two  combined  effects  apparently  favor  droplet  gasification  at 
high  pressure. 

To  summarize  the  overall  effect  of  ambient  flow  conditions  on  droplet  lifetime,  a  simple 
correlation  shown  in  Fig.  13  is  developed  as 


where  Re  is  the  Re3molds  number  based  on  the  initial  droplet  diameter  and  pr,Oa  the 
reduced  pressure  in  reference  to  the  critical  pressure  of  oxygen.  Equation  5.1  bears  a 
resemblance  to  the  classical  Ranz-Marshall  correlation  (1952)  for  droplet  vaporization  in 
convective  environments,  which  takes  the  form 


1  +  0.276ReV2pri/3 

The  Ranz-Marshall  correlation  is  only  applicable  to  low-pressure  conditions,  and  shows 
a  weaker  Reynolds-number  dependency. 

Because  of  lack  of  distinct  gas/liquid  interface  at  supercritical  conditions,  droplet  mo¬ 
tion  is  best  characterized  based  on  the  locus  of  the  center  of  gravity  which  can  be  obtained 
&om  the  first  moment  of  inertia  of  the  mass  confined  by  the  critical  mixing-temperature 
surface  in  reference  to  a  fixed  axis.  Figures  14(a)  and  (b)  present  the  results  for  various 
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ambient  freestream  velocities  and  pressures,  respectively.  The  droplet  travel  distance 

depends  on  both  the  net  momentum  transfer  and  the  time  duration  over  which  the  force 
acts  on  the  droplet.  For  instance,  at  p  =  100  afm,  although  the  convective  Reynolds 
number  aX  Uoo  =5  m/s  is  larger  than  that  at  1.5  m/s,  the  droplet  moves  farther  in  the 
low-velocity  case,  simply  due  to  its  longer  lifetime. 

The  motion  of  a  droplet  is  essentially  a  result  of  net  momentum  transfer  to  it.  Following 
common  practice,  the  drag  coefficient  is  defined  as  follows. 

Cd  =  1 - ^ -  (5.3) 

2pBoiUco  —  Uilrop)^vR^ 

where  {Uoo  —  Udrop)  is  the  droplet  relative  velocity,  and  Rc  the  radius  of  the  mass- 
equivalent  sphere  with  a  uniform  distribution  of  density  at  its  initial  value.  The  denomi¬ 
nator  in  Eq.  (5.3)  represents  the  aerodynamic  loading.  The  total  drag  acting  on  a  droplet 
includes  form,  &iction,  and  thrust  drag,  of  which  form  drag  is  usually  found  to  be  the 
dominant  component  (Ha3rwood  et  aL  1989).  Thrust  drag  arising  from  the  momentum 
transfer  with  the  ambient  flow  due  to  evaporation  is  always  negligibly  small.  Figure  15 
presents  the  time  histories  of  the  total  drag  coeflicient,  drag  force,  and  aerodynamic  load- 
P  =  100  atm  and  Uoo  =  5m/s,  where  the  asterisk  denotes  the  quantity  normalized 
by  the  initial  aerod3mamic  loading.  The  early  increase  in  is  caused  by  the  increasing 
form  drag  associated  with  droplet  deformation.  This  effect  is  then  counterbalanced  by 
the  strong  surface  blowing  due  to  gasification  which  weakens  both  friction  and  form  drag. 
The  total  drag  force  diminishes  to  zero  at  the  end  of  the  droplet  lifetime. 

The  drag  coefficient  of  an  evaporating  droplet  is  usually  smaller  than  that  of  a  hard 
sphere  at  the  same  Re3molds  number  (Chen  k  Yuen,  1976).  This  phenomenon  may  be 
attributed  to  the  modification  of  the  flowfield  and  transport  properties  in  the  vicinity  of 
the  droplet  surface.  The  gasification  process  often  leads  to  a  thickened  boundary  layer. 
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thereby  reducing  shear  stress  and  resultant  drag.  Several  researchers  (e.g.,  Renksizbulut 

&  Ha3rwood  1988;  Haywood  et  ai  1989;  Chiang  et  ai  1992)  niunerically  studied  this 
issue  by  solving  the  Navier-Stokes  equations  for  a  spherical  droplet  at  low  to  moderate 
pressures.  No  deformation  was  considered  to  simplify  the  analysis.  The  results  lead  to 
the  following  correlation 


Cd  = 


Cl 

(1  +  B)^ 


(5.4) 


where  Cl  is  the  reference  drag  coeflScient  for  a  hard  sphere(Putnam  1961)  and  takes  the 
form 


(5.5) 


The  coefficient  b  has  a  value  of  0.2  in  Renksizbulut ’s  model  (1988)  and  0.32  in  Chiang’s 
model  (1992).  The  Spalding  transfer  number  defined  below  is  adopted  to  account  for  the 
effect  of  evaporation  on  drag  coefficient. 


D  _  Cp{Too  ~  Ts) 

— 

Apparently,  the  above  correlation  becomes  invalid  in  the  present  study  of  supercritical 
droplet  gasification.  The  enthalpy  of  vaporization  Ahy  diminishes  to  zero  when  the 
thermodynamic  critical  state  is  reached,  rendering  a  singular  value  of  the  transfer  number. 

To  establish  a  generalized  correlation  for  supercritical  conditions,  the  instantaneous 
drag  force  was  calculated  throughout  the  entire  droplet  lifetime  for  ail  the  scenarios 
under  consideration.  The  data  was  then  correlated  through  the  use  of  a  transfer  number 
Bo  taking  into  account  the  rapid  transient  during  the  gasification  process,  defined  as 
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(5-') 

where  Ti  is  the  instantaneous  average  temperature  of  droplet,  and  Tc  the  critical  Tr.;-riT.g 
temperature.  Since  Bd  diverges  at  =  Ti  at  the  end  of  droplet  lifetime,  the  calcuUtion 
of  drag  force  was  terminated  when  %  -  Tt  becomes  less  than  1  K,  at  which  the  droplet 
residual  mass  is  usually  less  than  one  thousandth  of  the  initial  mass.  The  influence  on  the 
accuracy  of  data  reduction  is  quite  limited.  Following  the  procedure  to  Eq.  (5.4), 

a  correlation  shown  in  Fig.  16  for  LOX  droplet  drag  coefficient  U  obtained. 

=  {ifXy  (5-«) 

wh«e  a  ud  i  «e  sdected  lo  be  0.06  ud  1.592  lew^tiwly.  Tbe  d<te 

clusters  along  the  classical  drag  curve,  Eq.  (5.5),  in  the  low  Reynolds-number  region, 
but  deviates  considerably  at  high  Reynolds  numbers  (i.e.,  Re  >  10).  Although  a  shape 
fector  may  be  employed  to  account  for  this  phenomenon  arising  from  the  increased  form 
drag  due  to  droplet  deformation,  the  difficulty  of  calculating  this  factor  and  conducting 
the  associated  data  analysis  precludes  its  use  in  correlating  the  drag  coefficient 
Instead,  a  simple  correction  factor  Re°  ®  is  incorporated  into  Eq.  (5.8)  to  provide  the 
compensation.  The  final  result  shown  in  Fig.  17  is  given  below. 

QO  /JgO.3 

~  (l  +  aBT)i  S92  0  <  Re  <  300  (5.9) 

6.  Conclusion 

Gasification  of  LOX  droplets  in  supercritical  hydrogen  streams  has  been  analyzed 
numerically.  The  model  treats  the  complete  conservation  equations  in  an  axisymmetric 
coordinate,  and  accommodates  thermodynamic  nonideality  and  transport  anomaly.  The 
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analysis  enables  a  comprehensive  investigation  into  droplet  behavior  over  a  broad  range 

of  thermodynamic  state.  Various  key  processes  involved  in  droplet  gasification  in  a 
supercritical,  forced-convective  environment  were  identified,  including  flow  structures, 
droplet  d3mamics,  and  interphase  transport  phenomena.  In  addition,  a  parametric  study 
of  droplet  mass  and  momentum  transfer  as  a  function  of  ambient  pressure  (100-400 
atm)  and  Reynolds  number  (0-300)  was  conducted.  The  resultant  correlations  of  droplet 
lifetime  and  drag  coefficient  can  be  effectively  used  as  a  physical  submodel  in  high- 
pressure  spray  combustion  analyses. 
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Figure  1 .  Schematic  Diagram  of  Droplet  Vaporization  in  Forced  Convective  Environment. 
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Figure  2.  Comparison  of  Oxygen  Density  Predicted  by  the  BWR  Equation  of  State  and 

Measured  by  Sychev  et  ah  (1987). 
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Figure  3.  Rdative  Errors  of  Density  Predictions  by  Three  Different  Equations  of  State, 
Experimental  Data  from  Sychev  et  al.  (1987). 
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Figure  4.  Droplet  Reynolds  Number  as  a  Fnnction  of  Freestream  Velocity  and  Ambient 
Pressnre;  LOXIH2  System  with  Do  =  lOOpm  and  Too  =  lOOOiT. 
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Figure  5.  LOX  Droplet  Gasiftcation  in  Supercritical  Hydrogen  Flow;  p  =  100  atm,  Uoo  —  2.5  m/, 
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Figure  6.  LOX  Droplet  Gasification  in  Supercritical  Hydrogen  Flow;  p  =  400  atm,  Uoo  =  2.5  m/ 
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Figure  7.  LOX  Droplet  Gasification  in  Supercritical  Hydrogen  Flow;  p  =  100  atm,  Uoo  =  15  m/, 
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Figure  8.  LOX  Droplet  Gasification  in  Supercritical  Hydrogen  Flow  at  p  =  100  atm:  (a)  Spherical  Mode;  Uoo  =  0.2  m/s,  t=0.61  ms. 

(b)  Deformation  Mode;  Uoo  =  1.5  m/s,  t=0.61  ms.  (c)  Stripping  Mode;  Uoo  =  5  m/s,  t=0.17  ms,  (d)  Breakup  Mode;  Uoo  =  15  m/s,  t=0.17  ms. 
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Figure  9.  Effect  of  Ambient  Velocity  on  Evolution  of  Droplet  Critical  Surfaces,  LOX/H2  System  at  p  =  100  atm. 
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Figure  10.  Effect  of  Ambient  Pressure  on  Evolution  of  Droplet  Critical  Surfaces,  LOXIH2  System  at  Uoo  =  2.5  m/. 
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Figure  11.  Effect  of  Ambient  Velocity  on  Time  Variation  of  Droplet  Residual  Mass;  LOX/H2 

System  at  p  =  100  atm. 
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Figu^  12.  Effect  of  Ambient  Pressure  on  Time  Variation  of  Droplet  Residual  Mass; 
LOX/H2  System  at  Uoo  =  2.5  m/s. 
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Figure  13.  Correlation  of  Droplet  Lifetime  as  a  Function  of  Free-Stream  Reynolds  Number 

and  Pressure. 
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Figure  14.  Locus  of  Droplet  Center  of  Gravity,  LOX/H2  System:  (a)  Effect  of  Convective 
Velocity  at  p  =  100  atm.  (b)  Effect  of  Ambient  Pressure  at  Uoo  =  2.5  m/s. 
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Figure  15.  Time  Variation  of  Total  Drag  Force,  Aerodynamic  Loading,  and  Drag  Coefficient; 
LOX/H2  System  at  p=100  atm  and  Uoo  =  5  m/s. 
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This  paper  presents  a  comprehensive  theoretical  analysis  conducted  to  study  liquid  oxy¬ 
gen  (LOX)  droplet  vaporization  for  a  broad  range  of  ambient  conditions  in  either  hydro¬ 
gen  or  water /hydrogen  quiescent  environments.  An  important  effort  has  been  devoted 
to  the  estimation  of  trcinsport  and  thermodynamic  properties.  In  addition,  influence  of 
cross-diffusion  effects  (Dufour  and  Soret  effects)  is  studied,  showing  that  these  effects 
may  be  neglected  in  droplet  vaporization  studies  at  high  pressures.  Condensation  of 
gaseous  water  in  the  low-temperature  zone  close  to  oxygen  droplet  surface  is  taken  into 
accoimt  by  means  of  a  detailed  sub-model.  Predictions  show  that  condensation  occurs  in 
the  neighborhood  of  the  oxygen  droplet.  However,  the  condensed  particles  do  not  reach 
droplet  surface  but  axt  blown  away  by  the  convection  flow  induced  by  the  vaporization 
process.  Thus,  water  condensation  hcis  only  a  small  influence  on  the  vaporization  process. 
A  correlation  of  the  droplet  lifetime  valid  for  the  oxygen/hydrogen/ water  system  and  for 
ambient  conditions  such  that  supercritical  vaporization  regime  occurs  in  the  early  stages 
of  the  droplet  lifetime  is  presented. 


NOMENCLATURE 

a  Equation  of  state  constant  in  Eq.(2.8) 

A  Surface  area 

6  Equation  of  state  constant  in  Eq.(2.8) 

Spalding  number 
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Cp  Constant-pressure  specific  heat 

D  Droplet  diameter 

Dij  Binary  mass  diffusion  coefficient 

Dim  Effective  mass  diffusion  coefficient 

Dxi  Thermal  diffusion  coefficient  of  species  i 

e  Specific  total  internal  energy 

Ft  Thermophoretic  force 

Fv  Viscous  force 

G  Gibbs  free  energy  per  mole 

g  Gibbs  free  energy  per  unit  mass 

hi  Specific  enthalpy  of  species  i 

J  Homogeneous  nucleation  rate 

Kvap  Vaporization  kinetic  coefficient  in  Eq.(3.3) 

m  Droplet  mass  evaporation  rate 

N  Number  of  species 

p  Pressure 

Energy  diffusion  flux 
qi  Mass  diffusion  flux  of  species  i 

r  Radied  coordinate 

R  Droplet  radius 

Ru  Universal  gas  constant 

S  Saturation  ratio 

s  Specific  total  entropy 

V  Velocity 

Vs  Control  surface  moving  velocity 

V  Total  volume 

Vi  Diffusion  velocity  of  species  i 

Xi  Mole  fraction  of  species  i 

Yi  Mass  fraction  of  species  i 

W  Molecular  weight 
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Z  Compressibility  factor 

Greek  Symbols 

a  Thermctl  diffusivity 

a(T)  Soave  temperature  function  in  Eq.(2.8) 
A  Thermal  conductivity 

A  Molecular  mean  free  path 

//  Viscosity 

fit  Chemical  potential  of  species  i 

p  Density 

a  Surface  tension 

r  Droplet  lifetime 

(jj  Acentric  factor 

Subscripts 

c  Critic^ll  condition 

i  Species  i 

p  Condensed  particle 

r  Reduced  thermodynamic  property 

rf  Reference  fluid 

sat  Saturation  value 

oo  Ambient  condition 

Superscripts 
g  Gaseous  phase 

/  Liquid  phase 

0  Dilute  gas  limit 

Q  Time  derivative 

♦  Dimensionless  quantity 
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1.  Introduction 

Droplet  vaporization  is  a  physical  process  occurring  in  many  liquid-fueled  combustion 
devices.  Atomization  of  liquid  fuel  leads  to  the  formation  of  a  spray  of  droplets  that 
then  undergoes  a  sequence  of  vaporization,  ignition,  and  combustion  processes.  In  the 
past  few  decades,  combustor  performance  has  been  substantially  enhanced  by  increcising 
the  operating  temperature  and  pressure  in  the  thermodynamic  supercritical  state  of  the 
fuel  species.  The  conventional  models  developed  for  low-pressure  applications  could  no 
longer  been  applied.  Recently,  several  theoretical  works  have  been  devoted  to  the  under¬ 
standing  of  droplet  vaporization  and  combustion  under  high  pressure  conditions.  Both 
hydrocarbon  droplets  in  air,  see  for  example  Hsieh  et  al  (1991),  Shuen  ei  ah  (1992), 
Curtis  &  Parrel  (1992),  and  Jia  &  Gogos  (1993),  md  liquid  oxygen  droplets  in  hydrogen 
environments,  see  for  excimple  Litchford  &  Jeng  (1990),  Delplanque  &  Sirignano  (1991), 
Lin  et  al  (1994),  and  Daou  ei  al  (1995),  are  treated  comprehensively,  with  emphasis 
placed  on  the  effects  of  transient  diffusion,  interfacial  thermodynamics  with  dissolution 
of  gaiseous  components  in  the  liquid  phase.  The  work  of  Lin  ei  al  (1994)  appears  to  be 
the  most  comprehensive  to  date  because  of  the  attention  devoted  to  transport  properties 
modelization. 

The  present  work  is  concerned  with  the  vaporization  of  liquid  oxygen  in  either  pure 
hydrogen  or  mixed  hydrogen/water  environments  due  to  its  broad  applications  in  cryo¬ 
genic  rocket  engines  using  hydrogen  and  oxygen  as  propellants,  such  as  the  VULCAIN 
engine  (Ariane  5  main  stage  engine),  space  shuttle  main  engine  (SSME),  J-2  engine  (Sat¬ 
urn  V  3^^-stage  engine,  and  R?-10  engine  (?).  In  these  engines,  the  chamber  pressure 
and  temperature  overpass  the  critical  Vcilues  of  oxygen  and  hydrogen. 

In  the  present  study,  we  present  an  unified  treatment  of  transport  properties,  thermal 
conductivity  and  molecular  diflFusion  coefficient  valid  for  a  wide  range  of  thermodynamic 
states.  This  treatment  has  been  formulated  in  such  a  generalized  manner  that  the  results 
can  be  practically  applied  to  any  types  of  of  liquid-fueled  propellant  couples.  The  Soret 
and  Dufour’s  effects  have  been  estimated  carefully,  since  they  are  expected  to  be  signifi¬ 
cant  for  the  our  problem,  because  of  the  important  disparity  in  molecular  weight  between 
oxygen  and  hydrogen,  and  because  of  the  large  temperature  gradients  involved.  In  ad- 
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dition,  in  the  presented  analysis,  all  the  thermodynamic  properties  are  derived  from  the 
chosen  equation  of  state  in  a  very  self-consistent  manner.  This  self-consistent  treatment 
of  thermodynamic  properties  renders  the  numerical  scheme  very  efficient.  For  subcritical 
regimes  of  vaporization,  it  allows  the  use  of  a  kinetic-like  expression  for  evaluation  of  sub- 
critical  vaporization  rate,  which  avoids  any  iterative  matching  procedure.  All  the  previ¬ 
ous  features  constitute  prolongations  cmd  improvements  of  the  work  of  Lin  et  al.  (1994). 

Condensation  of  gaseous  water  diffusing  to  the  oxygen  droplet  surface  is  a  subject 
of  interest  for  many  years,  see  for  example  Litchford  &  Jeng  (1990)  and  Powell  (1988). 
This  phenomenon,  due  to  the  concomitant  effects  of  the  low  temperature  of  the  oxygen 
droplet  surface  and  of  the  diffusion  of  the  gaseous  water  towards  this  surface,  is  expected 
to  affect  the  surface  conditions  and  thus  the  vaporization  process-  All  the  numerical 
analyses  mentioned  previously  have  eluded  this,  and  we  can  only  make  reference  to  the 
simplified  analyses  of  Yoset  (1992),  Lafon  &  Prud’homme  (1994),  and  Lin  et  al  (1994). 
In  this  study,  a  comprehensive  analysis  of  the  water  condensation  process  is  carried  out 
by  estimating  all  the  different  mechcinisms  involved  in  this  process. 

Finally,  we  have  kept  in  mind  that  the  objective  of  any  single-droplet  vaporization 
study  is  to  provide  sub-models  used  in  spray  computations,  therefore  we  have  conducted 
a  large  series  of  computations  corresponding  to  a  variety  of  initial  and  ambient  condi¬ 
tions.  A  correlation  of  the  droplet  lifetime  has  been  deduced,  it  is  valid  for  the  oxy¬ 
gen/hydrogen/water  system  when  the  supercritical  vaporization  regime  occurs  in  the 
early  stage  of  the  droplet  lifetime. 

In  this  paper,  we  first  describe  the  theoretical  formulation  with  emphasis  on  transport 
and  thermodynamic  properties,  the  numerical  methods  used  are  thereafter  briefly  pre¬ 
sented.  The  model  is  first  applied  to  the  vaporization  of  oxygen  droplets  in  pure  hydrogen 
environments  and,  after  presenting  the  water  condensation  sub-model,  it  is  then  applied 
to  the  vaporization  of  oxygen  droplets  in  mixed  hydrogen/ water  environments. 

2.  Theoretical  Formulation 

We  consider  the  situation  where  a  droplet  initially  with  a  uniform  distribution  of  sub- 
criticail  temperature  is  suddenly  injected  in  a  supercritical,  quiescent  environment.  The 
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problem  is  entirely  governed  by  a  complete  set  of  conservation  equations.  Because  buoy¬ 
ancy  effects  and  forced  convection  are  not  considered,  the  problem  is  one-dimensional. 


2.1.  Conservation  Equations 

Badance  equations  are  valid  for  both  phases.  Let  us  consider  an  arbitrary  control  volume 
Va{t)  delimited  by  a  surface  Aa(t)  moving  in  an  absolute  referential  at  speed  Vs.  The 
mass  balance  equation  may  be  written  as  follows 

pdV  +  f  p(v  -  Vs)  •  dA  =  0  (2.1) 

at  JAa(t) 

The  momentum  and  energy  conservation  equations  are  simplified  with  the  following 
hypotheses.  The  velocities  involved  are  very  small,  so  that  total  energy  is  almost  equal 
to  interned  energy  and  viscous  dissipation  Cem  be  neglected.  The  momentum  equation  is 
then  reduced  to  a  space-constant  pressure  equation 


Vp  =  0 


(2.2) 


From  the  first  law  of  the  thermodynamics  applied  to  the  control  volume,  we  derive  the 
conservation  of  energy 


^  f  pe  dV  +  f  pe{v  —  •  dA  =  —  /  pv  '  dA—  [  q^-  dA 

Species  conservation  does  not  involve  any  simplifying  hypotheses 

f  pYidV+  f  pYi{v  -Vs)  •  dA  =  -  f  qj  •  dA 

JK(t)  JA,(t)  JA.(t) 


£ 

dt 


(2.3) 


(2.4) 


2.2.  Diffusion  Fluxes 

The  energy  and  mass  diffusion  fluxes,  qe  and  g,-,  in  Eqs.(2.3)  and  (2.4)  need  to  be  modeled. 
To  set  the  amalysis  in  its  most  general  form,  the  Soret  (thermal-diffusion)  effect  which 
accounts  for  heat  diffusion  due  to  concentration  gradients  is  considered  along  with  its 
reciprocal  Dufour  phenomenon.  The  Soret  effect  is  evidenced  by  the  observation  that 
light  molecules  tend  to  diffuse  towards  high-temperature  regions  and  vice  versa.  This 
phenomenon  may  be  significant  in  the  present  study  because  of  the  disparity  of  molecular 
weight  between  hydrogen  and  oxygen  and  the  steep  temperature  gradient  involved.  The 
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overall  diffusion  fluxes  take  the  form 

N  N  N 

q,  =  -Avr -'Zqi'hi+R.T^Y^  {Vi  -  Vj)  (2.5) 

vr 

qi  =  pYiVi  =  --pDimVYi  ^  Dn—  (2.6) 

where  Fourier’s  and  Pick’s  laws  are  used  to  approximate  the  ordinary  energy  and  mass 
diffusion,  respectively.  As  it  was  mentioned  in  Bird,  Steward  &  Lightfood  (1960),  for 
systems  involving  more  than  2  species,  Pick’s  law  is  only  an  approximation  of  the  molec¬ 
ular  diffusion  velocity  system  which  necessitates  the  use  of  a  mass  averaged  diffusion 
coefficient 

N 

Dim  =  (1  -  Xi)/  Y,{Xj/Dij)  (2.7) 

i^3 

2.3.  Property  Evaluation 

Transport  and  thermodynamic  properties  play  a  decisive  role  in  droplet  vaporization  pro¬ 
cess.  The  former  rule  mass  and  energy  transfer,  whereas  the  latter  affect  the  equilibrium 
conditions  and  energy  needed  for  phase  change.  Both  properties  dictate  the  conditions  of 
occurrence  of  supercritical  regimes.  Thus,  if  accurate  predictions  of  the  droplet  lifetime 
are  wanted,  precise  modelization  of  thermodynamic  behavior  and  transport  properties 
is  needed.  In  addition,  since  the  two  initial  phases  may  collapse  to  a  single  phase,  the 
property  evaluation  must  be  unified  over  the  entire  range  of  temperatures,  pressures  sind 
compositions. 

2.3.1.  Thermodynamic  Properties 

All  the  thermodynamic  properties  are  derived  from  the  Soave-Redlich-Kwong  equa¬ 
tion  of  state,  see  Graboski  k  Daubert  (1979).  It  represents  a  good  compromise  between 
computation  complexity  and  physical  accuracy,  as  was  mentioned  by  Lin  et  al  (1994). 
Pressure  may  be  expressed  in  terms  of  density  and  temperature 

_  pRuT  aoc  p^ 

{W  -hp)"^  W  {W  ^hp) 


(2.8) 
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where  a  and  b  accounts  for  attractive  forces  and  repulsive  forces  between  molecules, 

respectively.  Mixing  rules  to  compute  a  and  b  are  based  on  corresponding  state  principle 
N  N  N 

aa  =  ^^XiXj{l  -  kij)y/dJaJo^  ,  b  =  ^Xibi  (2.9) 

*  ;  i 

where  af,  bi  and  a,(T)  are  related  to  the  critical  properties  cind  acentric  factor  of  con¬ 


stituent  species  as  follows 


0.42748P2Tc? 


QM664RuTci 


(2.10) 


a,(T)  =  [l  +  fui  (l  -  x/r/Tc,)]  ^  with  fui  =  0.48  +  1.574w,-  -  0.176u;? 

All  the  thermodynamic  properties  required  in  the  analysis  come  from  derivation  of 
Eq.(2.8)  and  (2.9).  For  instance,  the  internal  energy  is  derived  from  the  following  relation 


(2.11) 


The  previous  integral  may  be  computed  to  give  the  following  analytic  formula 

e{p,  T,  Yi)  =  e°{T,  Yi)  +  ^  Log(l  +  bp/W)  (2.12) 

This  expression  of  internal  energy  may  be  written  as  an  explicit  function  of  temperature 
and  constituent  partial  density.  Then,  some  thermodynamic  considerations,  described  in 
Lafon  (1994),  enable  us  to  compute  easily  specific  enthalpy  of  each  species  using  direct 
derivation  the  Eq.(2.12) 


7  fdmh\  (dpe\  f213) 


To  compute  egression  of  chemical  potential,  first  expression  of  entropy  is  required. 
By  integration  of  the  relation 

)  =  *«(,, T,y)  +  ^io8Z+^'  P-w) 


an  analytical  expression  of  entropy  is  obtained 


+  iLog(l-V/«')  +  5^[^l  Log(l  +  v/»0 

L  J  Yi 


Supercritical  vaporization  of  liquid  oxygen  droplets 


9 


Chemical  potential  is  then  easily  evaluated  by  means  of 


(2.16) 


2.3.2.  Transport  Properties 

Transport  properties  determine  the  mass  and  energy  transport  which  dictate  the  in¬ 
terface  conditions  and  consequently  the  vaporization  rate  and  the  occurrence  of  critic^tl 
conditions.  The  major  properties  of  interest  are  thermal  conductivity,  mass  diffusivity 
and  thermal-diffusion  coefficient.  Classical  gas  kinetic  theory  fails  to  correctly  predict 
transport  properties  at  high  pressures.  However,  it  is  always  possible  to  estimate  in  the 
dilute-gas  limit  and  then  take  into  account  the  dense-fluid  effects.  Therm2d  conductivity 
is  generally  divided  into  two  contributions, 


Xm{p,T,Yi)  =  +  X‘Up>T,Yi) 


(2.17) 


A|,i(r,yi)  and  X*^{p,T,Yi)  represent  the  internal  contribution  and  the  translational  and 
collisional  contribution  to  the  thermal  conductivity,  respectively.  For  a  mixture  the  first 
term  may  be  evaluated  by  means  of  an  empirical  mixing  rule 

in  conjunction  with  the  Eucken’s  correlation  modified  for  polyatomic  gases 

A:.  =  1.32,i?(c°,-|^)  (2.19) 

where  the  low  pressure  estimations  of  viscosity  and  specific  heat  are  obtained  from  poly¬ 
nomial  adjustments  of  existing  experimental  data,  see  for  example  Reid,  Prausnitz  k 
Poling  (1987). 

The  collisional  and  translational  contribution  T,  Yi)  may  be  evaluated  using  the 
corresponding-state  method  proposed  by  Ely  &  Hanley  (1981)  and  (1983).  The  method 
is  able  to  predict  transport  properties  over  a  Icirge  domain  of  temperatures  and  pressures, 
from  compressed  liquid  to  dilute  gas.  The  scheme  assumes  that  the  configurational  prop¬ 
erties  of  a  single-phase  mixture  can  be  equated  to  those  of  a  hypothetical  pure  fluid,  which 
are  then  evaluated  via  corresponding  states  with  respect  to  a  given  reference  fluid.  It 
requires  the  characteristic  values  of  each  species:  critical  coordinates  and  acentric  factor. 
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Details  of  the  procedure  may  be  found  in  Ely  &  Hanley  (1983)  or  Hsiao  ei  al  (1995),  but 
let  us  briefly  present  the  guidelines  of  the  method  used  to  compute  thermal  conductivity. 

T,  Yi)  =  {Prf ,  Trf)  TxXx  (2.20) 

T\  and  X\  are  the  scaling  factor  eind  a  correction  factor  respectively.  The  latter  com¬ 
pensates  the  non-correspondent  behavior  of  the  mixture.  For  the  reference  fluid,  the 
collisional  and  translational  contribution  of  thermal  conductivity  is  divided  into  three 
parts 

^r/(Pr/,Tr/)  =  A°y(Tr/)  -f  AAearc(pr/ ,  ^r/)  +  AAcr*t(Pr/ ,  Tr/ )  (2.21) 

A°^(Tr/)  is  the  value  at  the  dilute  gas  limit  which  may  be  evaluated  by 

AAcxc(Pr/>rr/)  and  AAcr»t(pr/)2^r/)  Correspond  to  the  dense  fluid  correction  and  to  the 
critical  enhancement,  respectively.  They  may  be  computed  from  analytical  expressions 
fitted  to  experimental  data.  The  critical  enhancement  part  is  not  included  in  the  present 
analysis  for  two  reasons;  first,  properties  which  diverge  strongly  for  pure  components 
diverge  only  weakly  for  mixtures  as  discussed  by  Bruno  &  Ely  (1991),  second,  the  crit- 
ical  enhancement  of  the  thermal  conductivity  has  not  been  observed  experimentally  for 
mixtures. 

Let  us  now  present  the  estimations  techniques  used  for  the  molecular  diffusion  coef¬ 
ficient.  For  gaseous  phase,  the  approach  proceeds  in  two  steps.  A  low  pressure  binary 
mass  diffusivity  is  obtained  using  the  Chapman-Enskog  theory  in  conjunction  with  the 
Lennard-Jones  intermolecular  potential-energy  function,  see  for  example  Reid,  Prausnitz 
&  Poling  (1987).  The  coefficient  is  corrected  with  the  correlation  of  due  to  the  work  of 
Takahashi  (1974)  to  take  into  account  high  pressure  effects.  However  this  method  fails 
to  predict  liquid  molecular  diffusion  coefficient  which  is  best  estimated  with  Wilke-and- 
Chang  procedure.  As  supercritical  regime  is  reached  liquid  estimation  technique  is  still 
applied  for  temperatures  below  the  mixture  critical  temperature.  These  methods  ^lre 
matched  at  the  interface  by  using  a  correcting  factor  for  the  liquid  estimation. 

Concerning  the  thermal  diffusion  coefficients,  much  less  material  is  available  in  the 
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literature,  especially  about  the  pressure  effects  on  these  coefficients.  Pressure  effects 
are  expected  to  be  of  the  same  order  of  magnitude  as  on  binary  molecular  diffusion 
coefficients.  However,  to  estimate  these  coefficients,  we  use  the  theoretical  developments 
of  Dixon-Lewis  (1968).  This  work  is  based  on  extension  of  the  Chapman-Enskog  kinetic 
to  polyatomic  gases.  Since  this  method  is  valid  only  in  the  low  pressure  region,  over¬ 
estimations  of  the  thermal  diffusion  coefficients  are  expected. 


2.4.  Interfacial  Boundary  Conditions 

Before  the  surface  temperature  reaches  the  mixture  critical  value,  thermodynamic  phase 
equilibrium,  characterized  by  a  minimum  of  the  Gibbs  function,  is  assumed  to  prevail  at 
the  interface.  If  surface  tension  is  neglected,  the  equilibrium  relations  may  be  written  as 


j'g 

=  r' 

pO 

=  p' 

= 

(2.22) 

A 

The  conservation  of  mass,  species  and  energy  across  the  interface  are 


m  =  p(v  ~  R)i4|  =  p(v  —  K)A 

\r=R+ 

r=R- 

(2.23) 

ihi  =  [ml^-  +  =  [mi;-  + 

(2.24) 

N 

(2.25) 

where  iZ_  and  denotes  the  conditions  at  the  interface  on  the  liquid  side  and  gaseous 
side,  respectively. 


3,  Numerical  Methods 

The  set  of  equations,  Eqs.  (2.1)-(2.4),  is  solved  using  a  fully  implicit  scheme.  Bal¬ 
ance  equations  are  discretized  on  a  time  varying  grid  with  a  finite  volume  formulation. 
Two  different  scenarios  occur  depending  upon  whether  the  vaporization  regime  is  sub- 
critical  or  supercritical.  However,  in  both  cases  the  grid  follows  the  time  regression  of 
droplet  interface.  For  supercritical  regimes  the  droplet  interface  is  either  defined  by 
critical  isotherm  or  critical  composition  line.  To  perform  a  fully  implicit  treatment  of 
the  problem  using  time- varying  grid  cell,  the  volume  of  each  grid  cell  is  included  in  the 
quantities  computed  at  each  time  step  by  inversion  of  the  linear  system  introduced  by 
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the  implicit  formulation.  The  variation  of  a  grid  cell  is  determined  by  the  difference  of 
its  limiting  surface  velocities.  If  the  location  of  each  grid  node  is  only  a  function  of  the 
instantaneous  droplet  radius,  then  the  volume  variation  of  each  grid  cell  may  be  easily 
linked  to  the  volume  variations  of  the  adjacent  cells,  same  thing  for  the  surface  veloci¬ 
ties.  Droplet  surface  velocity  is  given  by  Eq.  (2.23)  for  subcritical  vaporization  regimes, 
or  by  cissuming  that  the  interface  corresponds  to  the  mixture  critical  temperature  line  or 
critical  composition  line  for  supercritical  vaporization  regimes.  Detsdls  of  the  numerical 
methods  are  presented  in  Lafon  (1994).  The  following  vector  summarize  the  quantities 
computed  at  each  time  step  for  each  grid  cell 

Q  =  {pV,flvA,peV,piV,Vf  (3.1) 


All  the  quantities  required  in  the  resolution  are  expressed  with  those  contained  in  the 
previous  vector  using  a  linearization  technique.  This  technique  gives  for  the  expression 
of  temperature 


2^+1  ^  j,  (gn+1  _  gn) 


(3.2) 


where  the  expression  of  the  vector  VqT  is  derived  from  thermodynamic  considerations 
involving  Eq.(2.8)  and  its  derivatives,  see  Lafon  (1994)  for  details.  Time  step  is  reduced 
when  temperature  variations  exceed  a  certain  value,  so  that  linearization  techniques  is 
still  applicable. 

The  treatment  of  interfacial  conditions  requires  a  special  procedure.  Liquid  and 
gaseous  conditions  at  the  interface  are  supposed  to  be  those  of  the  two  cells  adjacent 
to  the  surface.  Interface  thermal  conductivity  is  artificially  increased.  This  procedure 
does  not  alter  the  local  energy  flux  value  and  thus  reduced  the  local  temperature  gra¬ 
dient  so  that  temperature  difference  between  liquid  and  gaseous  pheise  vanishes  to  zero. 
Then  the  trctnsfer  of  each  species  involved  in  the  thermodynamic  equilibrium  between 
the  two  phases  is  realized  by  using  the  following  pseudo-phenomenological  or  kinetic-like 
relationship 

m  =  Kvap,i  {(i’i  -  /if)  (3.3) 


This  expression  is  comparable  to  phenomenological  laws  such  as  Fourier’s  law  for  heat 
conduction  or  Pick’s  law  for  mass  diffusion.  The  flux  is  proportional  to  the  difference 
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of  potential  via  the  diffusion  coefficient.  To  ensure  the  thermodyneimic  equilibrium,  the 
value  of  the  vaporization  kinetic  coefficient  is  set  so  that  the  difference  between  chemical 
potential  is  as  small  as  wanted.  This  technique  requires  axi  implicit  expression  which  is 
obtained  by  means  of  a  relation  similar  to  Eq.(3.2)  for  chemical  potential.  This  method 
does  not  necessitate  amy  sub-iterative  process,  but  it  requires  an  expression  of  chemical 
potential  derivatives  as  accurate  and  as  consistent  as  possible.  If  not,  computations 
will  diverge.  That  why  all  thermodynamic  properties  are  only  derived  from  derivations 
of  Eqs.(2.8)  and  (2.9).  Once  the  thermodynamic  equilibrium  is  initiated,  the  interface 
evolves  on  a  thermodynamic  equilibrium  trajectory.  This  evolution  is  controlled  by  heat 
amd  matss  fluxes  at  the  interface.  As  soon  ats  the  interface  reaches  the  mixture  critical 
conditions,  this  procedure  is  no  longer  applied,  and  the  anailysis  switches  to  a  single-phase 
approach. 

4.  Vaporization  in  Hydrogen 

The  model  is  first  applied  to  the  vaporization  of  liquid  oxygen  droplet  in  hydrogen 
environments  with  a  broad  range  of  initial  amd  ambient  conditions. 

4.1.  Binary  Equilibrium  and  Criticality 

Let  us  consider  the  binary  equilibrium  involved  in  this  pairticular  case  at  the  droplet 
interfaw:e.  The  variance  of  the  thermodynamic  equilibrium  is 

=  JV-f.2-(?^  (4.1) 

Only  two  species  (oxygen  and  hydrogen)  and  two  phases  (liquid  amd  gas)  are  present, 
so  that  the  vauriance  is  2.  Furthermore,  the  momentum  balance  over  the  domain,  amd 
especiadly  at  the  interface,  prescribes  a  constant  pressure.  As  a  consequence  the  reduced 
variance  is  unity.  Evolutions  of  composition  and  temperature  at  the  interface  are  linked 
by  the  thermodynamic  equilibrium.  Trajectories  of  the  equilibrium  are  reported  in  Fig.  1. 

For  subcritical  pressures  hydrogen  dissolution  in  the  liquid  phase  is  small  and  almost 
negligible.  As  the  pressure  increases  above  the  critical  pressure  of  pure  oxygen  (49.7  a^m), 
the  dissolution  effect  becomes  more  and  more  significant.  For  a  fixed  pressure,  dissolution 
is  maodmum  when  the  two  phases  reach  their  critical  conditions,  and  increases  generally 
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with  pressure.  Thus,  for  pressure  equal  to  300  atm,  critical  hydrogen  molar  composition 
is  close  to  0.5.  However,  since  there  is  a  big  disparity  in  molecular  weight  between  both 
species,  hydrogen  critical  mass  compositions  are  still  close  to  zero.  Nevertheless,  hydrogen 
properties  are  different  enough  from  the  oxygen  ones  to  affect  substantially  the  mixture 
properties.  As  an  example,  enthalpy  of  vaporization  is  reported  for  various  pressures  in 
Fig.  2.  Dissolution  makes  enthalpy  of  vaporization  a  function  of  both  temperature  and 
pressure. 

Let  us  discuss  now  about  an  important  feature,  which  concerns  the  properties  of  the 
mixture  critical  point.  For  a  pure  component  the  critical  point  is  the  locus  of  a  loss  of 
mechanical  stability;  the  first  and  second  derivatives  of  pressure  with  respect  to  den¬ 
sity  vanish  to  zero.  This  implies  a  strong  divergence  of  the  specific  heat  as  the  critical 
point  is  approached.  For  a  mixture,  the  criticality  criteria  is  somehow  different.  As 
its  critical  conditions  are  approached,  the  mixture  lose  its  material  stability,  a  phe¬ 
nomenon  which  may  be  characterized  for  a  binary  mixture  by  the  relation  as  mentioned 
by  Bruno  k  Ely  (1991) 


All  cubic  equations  of  state  use  mixture  rules  based  on  the  corresponding  state  prin¬ 
ciple.  The  mixture  properties  are  reported  to  the  properties  of  a  hypothetical  pure 
component  by  means  of  the  mixture  rules.  Thus,  for  the  mixture,  a  pseudo-critical  point 
exists;  it  has  all  the  characteristics  of  a  pure  component  critical  point,  and  is  therefore 
a  locus  of  a  strong  divergence  of  the  specific  heat.  However,  let  us  consider  a  binary 
mixture  which  pressure  is  above  the  critical  pressure  of  the  less  volatile  species  (oxygen), 
but,  which  temperature  is  below  the  mixture  critical  value.  When  temperature  progres¬ 
sively  increases,  it  will  reach  its  mixture  critical  point  before  the  pseudo-critical  point. 
Since  this  point  is  different  from  the  pseudo-critical  point  but  yet  not  too  far  from  it, 
only  a  weak  divergence  on  the  specific  heat  is  observed.  This  divergence  vanishes  when 
the  pressure  is  far  above  the  critical  pressure. 

When  interface  reaches  its  critical  conditions,  no  distinction  between  the  two  phases 
remains.  The  critical  interface  is  then  claissically  defined  as  the  critical  isotherm,  and 
it  is  still  the  locus  of  important  density  gradients.  Another  definition  based  on  critical 
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composition  may  be  considered.  This  choice  can  be  justified  by  the  fact  that  it  may  better 
represent  the  mixing  process  involved  in  supercritical  vaporization.  Application  of  one 
or  the  other  definition  will  lead  to  different  results  mainly  because  the  Lewis  number 
is  not  unity  (See  next  sections).  The  former  definition  seems  the  most  adequate  for 
comparison  with  experimental  data,  in  which  droplet  regression  is  tracked  by  means  of 
measurement  techniques  sensitive  to  temperature  or  density  gradients  -shadowgraphy  is 
used  by  Hsiang  &  Feath  (1993)-.  For  use  in  spray  analysis,  in  which  vaporization  process 
is  treated  as  a  source  term  of  mass  and  a  sink  term  of  energy  in  balance  equations,  the 
definition  of  the  critical  interface  based  on  temperature  seems  also  adequate.  The  mixing 
process  experienced  by  the  vaporized  species  is  then  hemdle  by  the  gas-mixing  model  of 
the  spray  model.  However,  in  the  following  study,  both  definition  based  on  temperature 
and  definition  based  on  composition  are  considered. 


4.2.  Discussion  of  Results 

A  series  of  calculations  have  been  carried  out  to  study  the  effects  of  eunbient  pressure 
and  temperature  on  droplet  behavior  in  both  the  subcritical  and  supercritical  vapor¬ 
ization  regimes.  In  this  pure  vaporization  case,  the  problem  is  entirely  controlled  by 
diffusion,  as  was  reported  by  Daou  et  al  (1995),  dimensional  analysis  shows  that  the 
initial  droplet  diameter  is  the  only  length  scale  and  that  the  droplet  lifetime  is  propor¬ 
tional  to  the  squared  initial  droplet  diameter.  Thus  only  droplets  with  an  initial  diameter 
jDo  =  100  are  considered  in  the  present  analysis. 

The  case  with  an  ambient  hydrogen  temperature  of  1000  K  is  first  considered  to  provide 
a  baseline  data  set.  Figure  3  shows  the  time  variations  of  droplet  surface  temperature 
at  various  pressures.  Three  different  scenarios  must  be  reported.  First,  for  low  pressures 
(p  =  10  atm),  the  surface  temperature  rises  suddenly  and  reaches  a  constant  value,  a  con¬ 
dition  referred  by  Hsieh  et  ah  (1991)  as  the  pseudo  wet-bulb  state.  The  constant  value  is 
slightly  lower  than  the  oxygen  boiling  temperature  because  of  the  presence  of  hydrogen 
on  the  gaseous  side  of  the  interface.  For  higher  pressures  (p  =  50  atm)  surface  temper¬ 
ature  rises  continuously.  The  pseudo  wet-bulb  state  disappears,  and  the  vaporization 
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process  becomes  transient  in  nature  during  the  entire  droplet  lifetime.  For  p  =  100  atm^ 
the  droplet  surface  even  reaches  its  critical  state  at  t  =  1  ms. 

In  Figs.  4,  5,  and  6,  the  instantaneous  distributions  of  the  mixture  specific  heat,  the 
mixture  thermal  diffusivity,  and  the  Lewis  number,  respectively,  are  reported  at  different 
times  for  p  =  100  afm,  Too  =  1000  K  and  To  =  90  K.  These  figures  show  the  extreme 
space  and  time  variations  of  the  transport  properties.  In  addition,  figure  4  clearly  shows 
the  weak  divergence  of  the  specific  heat  around  the  critical  conditions.  The  Lewis  number, 
defined  as  £e  =  X/pCpD  and  reported  in  Fig.  6,  characterizes  the  relative  significance  of 
thermal  diffusion  to  mass  diffusion. 

Figure  7  presents  the  time  evolution  of  the  reduced  droplet  diameter  squared  at  various 
pressures.  For  supercritical  regimes,  we  follow  the  critical  temperature  regression.  For 
low  pressures,  the  behavior  follows  that  predicted  by  the  qucisi-steady  theory  assumption. 
However,  at  high  pressures,  the  variation  becomes  nonlinear  with  respect  to  time.  The 
slope  dD^/dt  is  even  almost  infinite  at  the  beginning  of  the  process,  because  of  the  infinite 
temperature  gradient  at  the  droplet  surface. 

In  Fig.  9,  the  droplet  lifetime  is  reported  for  various  ambient  temperatures  as  a  func¬ 
tion  of  the  reduced  pressure  Pr  =  p/Pco^*  figure  use  logarithm  scales  on  both  axis. 
The  plots  exhibit  two  different  slopes,  depending  on  whether  the  vaporization  regime  is 
subcritical  or  supercritical.  For  subcritical  regimes,  the  droplet  lifetime  decreases  with 
ambient  pressure,  mainly  because  of  the  reduced  enthalpy  of  vaporization.  For  supercrit- 
icaJ  regimes,  the  droplet  behavior  is  then  dominated  by  heat  diffusion.  The  decrease  of 
critical  mixture  temperature  associated  with  the  increase  of  thermal  conductivity  with 
pressure  leads  to  a  decreasing  droplet  lifetime  with  pressure.  One  interesting  feature  is 
that  no  abrupt  change  is  noticeable  at  the  transition  despite  the  vanishing  of  enthalpy 
of  vaporization.  As  critical  conditions  are  approached,  for  pressures  slightly  above  the 
critical  pressure  of  pure  oxygen,  specific  heat  of  the  mixture  exhibits  a  divergence,  see 
Fig.  4.  The  closer  the  pressure  is  to  the  critical  pressure  of  pure  oxygen,  the  stronger  is 
this  divergence.  Then,  as  criticeil  conditions  are  approached,  although  the  heat  consumed 
by  vaporization  process  becomes  smaller,  the  heat  needed  to  incre3ise  the  droplet  surface 
to  its  critical  value  becomes  bigger.  These  two  conflicting  effects  tend  to  balance. 
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Figure  8  presents  the  time  evolution  of  the  reduced  droplet  diameter  squared  at  var¬ 
ious  pressures.  Contrary  to  the  previous  figure,  the  critical  interface  is  now  defined  as 
the  critical  composition  line.  This  lead  to  results  substantially  different  from  the  ones 
presented  in  Fig.  7.  In  Fig.  6,  one  can  notice  that  the  Lewis  number  is  smaller  than 
unity  in  the  gaseous  region  close  to  droplet  surface.  Thus,  heat  transfer  penetrates  faster 
the  droplet  interior  than  molecular  transfer,  all  the  more  so  that  pressure  is  high.  Then, 
critical  composition  line  is  blown  away  from  the  liquid  core.  Eventually  it  reaches  re¬ 
gions  where  molecular  diffusion  is  more  efficient,  higher  temperature  regions  leading  to  a 
Lewis  number  bigger  than  unity.  For  supercritical  regimes,  the  final  stage  of  the  droplet 
lifetime  is  characterized  by  a  high  regression  rate.  For  low  ambient  temperature,  where 
favorable  conditions  for  efficient  molecular  diffusion  are  more  difficult  to  find,  this  mech¬ 
anism  may  lead  to  an  increase  of  the  droplet  lifetime  for  high  pressure,  with  a  minimum 
occurring  for  an  intermediate  pressure  greater  than  oxygen  critical  pressure.  This  feature 
is  reported  in  Fig.  10,  which  present  the  droplet  lifetime  as  a  function  of  the  reduced  pres¬ 
sure  pr  =  pIPc02  various  ambient  temperatures.  With  the  critical  interface  based  on 
critical  composition,  no  significant  change  is  noticeable  at  the  transition  from  subcritical 
to  supercritical  regimes. 

Let  us  briefly  comment  on  the  influence  of  the  cross-diffusion  effects,  namely  the  Soret 
and  Dufour  effects.  Two  series  of  computations,  with  and  without  taking  into  account 
these  effects,  have  been  conducted.  These  computations  are  reported  in  Figs.  11  and  12. 
It  is  important  to  note  that  both  axis  use  linear  scale.  The  influence  of  thermal  diffusion 
is  noticeable  for  subcritical  pressures  as  well  as  for  supercritical  pressures  when  the 
critical  interface  is  defined  as  the  critical  composition  line.  The  most  visible  effect  of 
thermal  diffusion  is  the  fact  that  heavier  molecules  are  attracted  by  low  temperature 
zones  and  vice  versa  for  lighter  molecules.  In  our  case,  hydrogen  is  attracted  by  hotter 
regions  while  oxygen  tends  to  stay  in  cooler  regions.  For  the  subcritical  case,  since  the 
hydrogen  diffusion  flux  at  the  interface  is  smaller,  the  thermodynamic  equilibrium  is 
displaced  to  smaller  temperatures  (See  Fig.  1).  Figure  2  demonstrates  that  it  leads  to 
a  slightly  higher  vaporization  enthalpy.  This  phenomenon  tends  eventually  to  increase 
the  droplet  lifetime.  For  supercritical  regimes,  when  the  criticzil  interface  is  based  on 
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critical  composition,  the  rate  of  vaporization  is  controlled  by  the  species  transport.  The 
droplet  lifetime  is  affected  by  cross-diffusion  effects.  Once  again,  the  droplet  lifetime  is 
increased  in  comparison  to  the  case  without  thermal  diffusion.  The  Soret  effect  tends  to 
keep  oxygen  in  the  cold  region,  the  critical  composition  line  corresponds  then  to  a  lower 
temperature  where  species  transport  is  globally  less  efficient. 

Our  estimation  of  thermcd  diffusion  coefficients  should  be  restricted  to  low  pressures. 
Such  low-pressure  estimations  lead  to  thermal  coefficients  independent  of  pressure,  which 
is  more  likely  questionable.  In  fact,  the  dependency  of  thermal  diffusion  coefficients  with 
respect  to  pressure  must  be  roughly  the  same  that  for  molecular  diffusion  coefficients. 
Thus,  in  our  study,  effect  of  cross-diffusion  have  been  over-estimated  for  moderate  and 
high  pressures.  Furthermore,  since  this  effect  tends  to  increase  the  droplet  lifetime  by 
less  than  a  few  percents,  we  do  not  consider  it  in  the  following  presented  results. 

4.3.  Limits  of  Critical  Regimes 

Since, for  subcritical  regimes,  evolutions  of  temperature  and  compositions  are  linked  by 
the  thermodynamic  equilibrium,  the  evolution  of  interface  conditions  is  monitored  both 
by  heat  transfer  (evolution  of  temperature)  and  species  transport  (evolution  of  compo¬ 
sition).  The  limit  of  occurrence  of  supercritical  regimes,  in  terms  of  ambient  conditions, 
is  then  a  strong  function  of  the  estimation  of  transport  properties.  In  Fig.  13,  the  limit 
ambient  temperature,  at  which  critical  conditions  are  reached,  is  plotted  for  various  pres¬ 
sures,  and  two  initial  droplet  temperatures,  90  and  110  if,  respectively.  The  dashed  lines 
correspond  to  a  supercritical  regime  reached  after  1  %  of  droplet  mass  is  vaporized,  the 
solid  line  corresponds  to  a  supercritical  regime  reached  at  the  end  of  the  droplet  lifetime. 
For  a  pressure-temperature  couple  located  in  the  region  above  the  solid  line,  supercritical 
regime  occurs  during  the  droplet  lifetime.  For  pressures  above  the  double  of  the  crit¬ 
ical  pressure  of  pure  oxygen,  supercritical  regimes  are  reached  almost  instantaneously 
whatever  the  ambient  temperature  is.  There  are  two  major  explanations  of  this  phe¬ 
nomenon.  First,  as  pressure  increases,  the  mixture  critical  temperature  decreases  (See 
Fig.  1).  Consequently,  a  droplet  can  reach  its  critical  conditions  more  easily  at  a  higher 
pressure.  Second,  the  molecular  diffusion  coefficient  is  a  decreasing  function  of  pressure. 
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Then,  since  hydrogen  will  diffuse  less  at  higher  pressures  to  the  droplet  surface,  critical 
conditions  will  be  easier  to  reach.  There  is  another  interesting  application  of  the  last 
point:  if  the  hydrogen  ambient  composition  is  decreased  or  if  oxygen  ambient  composi¬ 
tion  is  increased,  by  considering  for  instance  a  mixture  of  hydrogen  and  oxygen,  the  limit 
temperature  for  supercritical  regime  will  decrease  with  respect  to  the  ones  presented  in 
Fig.  13.  To  conclude  with  this  point,  the  comparison  of  results  obtained  for  two  different 
initial  droplet  temperatures  demonstrated  that  the  limit  temperature  of  occurrence  of 
supercritical  regime  is  also  a  function  of  the  initial  droplet  temperature;  the  cooler  the 
droplet  is  initially,  the  hotter  ambient  gas  must  be  in  order  to  reach  supercritical  regime. 
This  effect  is  stronger  for  relatively  low  pressures. 

Umemura  (1986) 


developed  a  theoretical  model  of  droplet  vaporization  under  supercritical  conditions. 
The  model  was  one-dimensional  with  a  planar  geometry,  which  enabled  the  author  to  find 
a  self-similar  solution  to  the  transient  problem.  With  a  reasoning  involving  molecular 
and  heat  fluxes  at  the  interface  as  well  as  the  slope  of  the  equilibrium  curves  of  Fig.  1, 
Umemura  stated  that  any  thermodynaunic  model,  in  order  to  be  consistent,  requires  that 
molecular  diffusion  coefficient  vanishes  at  the  mixture  critical  point.  This  phenomenon, 
which  was  observed  experimentally  by  Sengers  (1972)  may  be  understood  simply  as  fol¬ 
lows.  For  sub  critical  regimes,  as  thermodynamic  equilibrium  exists,  a  difference  in  liquid 
and  gaseous  composition  at  the  interface  exists.  Since  interface  thickness  is  assumed  to 
be  infinitely  small,  the  gradient  of  composition  is  then  infinite.  As  molecular  flux  remains 
finite  one  may  consider  that  the  molecular  diffusion  coefficient  must  vanish  across  the  in¬ 
terface.  This  feature  must  remans  true  as  temperature  is  increaised  to  the  critical  value. 
However,  since  there  is  no  accurate  modelization  of  this  phenomenon  in  the  literature, 
we  did  not  implement  it  in  our  numerical  model.  From  the  previous  discussion,  we  can 
easily  understand  that  this  phenomenon  influences  the  limit  values  of  temperature  of  su¬ 
percritical  regime  occurrence  reported  in  Fig.  13.  These  latter  values  are  then  expected 
to  be  smaller  than  the  ones  presented. 
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4.4.  Approximate  Analysis 


In  this  section,  the  attention  is  focused  on  the  estimation  of  the  droplet  lifetime  re¬ 
stricted  to  supercritical  vaporization  regimes,  when  this  regime  occurs  soon  after  the 
liquid  droplet  have  been  exposed  to  the  hot  gaseous  ambient  atmosphere,  that  is  when 
ambient  pressure  eind  temperature  are  located  above  the  dashed  lines  of  Fig.  13.  The 
critical  interface  is  based  on  the  mixture  critical  temperature.  In  this  case,  as  mentioned 
in  the  previously,  the  problem  is  entirely  driven  by  heat  conduction.  In  order  to  obtain  an 
analytical  solution,  a  simplified  analysis  of  this  problem  may  be  conducted.  Simplifying 
hypotheses  are  thus  needed:  thermal  conductivity,  specific  heat  and  thermal  diffusivity 
a  =  X/pCp  are  constant  in  both  space  and  time.  In  addition,  the  Stefan  flow  induced  by 
vaporization  is  assumed  to  be  negligible.  Thus,  the  problem  reduces  to  the  resolution  of 
the  following  partial  derivative  equation 


JlA  (  2^^ 

dt  ^  dr  \  dr  ) 


(4.3) 


Using  thermal  diffusivity,  initial  droplet  radius,  initial  droplet  temperature,  and  ambient 
temperature,  we  can  substitute  in  this  equation  the  dimensionless  quantities:  r*  =  r/iZoj 


t*  =  ta/Rl,  and  =  (Too  -  T)/(Too  -  To),  to  obtain 


dt^  r*2  dr*  V  / 


(4.4) 


Considering  the  initial  temperature  distribution,  for  =  0 

rT*(r%0)  =  l  ,  forr"<l 

lT^(r*,0)  =  0  ,  forr*>l, 

The  solution  of  Eq.(4.4)  was  given  by  Rosner  (1967).  At  that  time  Rosner  wcis  interested 
in  the  burning  time  of  droplets  experiencing  supercritical  vaporization.  The  equation 
concerned  the  mass  fraction,  and  by  tracking  the  stoichiometric  location,  it  gave  pre¬ 
diction  of  the  burning  time.  Nevertheless,  the  solution  for  reduced  temperature  may  be 
given  by 


(4.5) 
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If  the  classical  definition  of  the  critical  interface  based  on  the  mixture  critical  temperature 
is  used,  from  the  previous  formulation,  the  reduced  droplet  lifetime  r*  is  demonstrated 
to  be  only  a  function  of  the  reduced  critical  temperature,  T*  =  (To©  —  Tc)/(Too  —  To), 
and  thus  solution  of  the  equation 

=  («) 

4*5.  Final  Droplet  Lifetime  Results 

A  large  series  of  numerical  simulations  have  been  conducted  for  various  ambient  tem¬ 
peratures  and  pressures  and  for  two  initial  droplet  temperatures:  To  =  90  and  110  JRT, 
respectively.  The  former  corresponds  to  the  lower  limit  of  LOX  injection  temperature 
in  many  operational  rocket  engines,  the  latter  is  used  to  determine  the  eflfect  of  initial 
droplet  temperature  and  is  slightly  under  critical  temperature  for  pressure  of  300  atm. 
Ambient  temperatures  were  chosen  in  the  range  500  <  Too  <  2500  AT,  and  ambient 
pressures  Pc  <P  <  300  atm.  Of  course  some  of  the  computations  may  not  meet  the  re¬ 
quirement  of  quick  criticality  and  will  present  significant  differences  from  the  other  ones. 
This  particular  case  is  however  representative  of  most  of  the  cases:  liquid  and  gaseous 
thermal  diffusivity  differ  by  several  orders  of  magnitude.  It  is  due  to  the  difference  be¬ 
tween  oxygen  and  hydrogen  behavior,  as  well  as  to  the  temperature  difference  between 
both  phases.  It  appears  clearly  that  the  characteristic  heat-up  process  in  the  liquid  phase 
is  the  limiting  process  of  droplet  vaporization.  Using  this  observation,  Figure  14  shows, 
for  all  the  performed  computations,  the  dimensionless  droplet  lifetime  r*  =  Tao/i?o  as 
a  function  of  the  reduced  critical  temperature  T*,  where  denotes  the  initial  liquid 
thermal  diffusivity.  The  solid  line  corresponds  to  theoretical  predictions  from  resolution 
of  Eq.(4.6).  The  filled  dots  correspond  to  computations  that  do  not  meet  the  fast  criti¬ 
cality  requirement,  ambient  pressure  and  temperature  located  under  the  dashed  lines  of 
Pig.  13. 

Unfortunately  no  line  meets  all  the  data  points,  clearly  because  it  involves  many  fea¬ 
tures  which  are  not  included  in  the  simplified  analysis.  First,  the  convection  induced 
by  the  vaporization  process  Weis  not  considered  in  the  analysis.  This  effect  reduced  the 
overall  heat  flux  arriving  at  the  droplet  surface  and,  if  it  wats  included  in  the  previous 
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theory,  it  would  shift  up  the  solid  line  of  Fig.  14.  However,  this  effect  is  more  likely  to 
be  small.  Since  the  process  is  entirely  driven  by  heat  diffusion,  the  space  distribution  of 
thermal  diffusivity  must  be  determincint,  all  the  more  that  liquid  and  ambient  thermal 
diffusivity  differ  by  several  orders  of  magnitude.  Figure  5  demonstrates  that  the  thermal 
diffusivity  at  the  droplet  surface  is  very  close  to  the  liquid  thermal  diffusivity;  the  same 
remark  may  be  done  for  thermal  conductivity.  The  heat  flux  arriving  at  the  droplet 
surface  is  the  product  of  the  interface  thermal  conductivity  multiplied  by  the  local  tem¬ 
perature  gradient.  The  interface  thermal  conductivity  is  almost  unchanged  during  the 
entire  droplet  lifetime,  but  the  interface  temperature  gradient  varies.  It  may  roughly  be 
estimated  as  follows 


VT| 


(4.7) 


As  the  ambient  thermal  diffusivity  increases,  the  interface  temperature  gradient  is  re¬ 
duced  so  that  the  droplet  lifetime  is  expected  to  be  longer  with  respect  to  the  case  where 
thermail  diffusivity  is  space  constant.  It  is  then  logical  to  correct  the  results  presented  in 
Fig.  14  by  a  factor  dependent  of  the  ratio  ttoo/ofo-  By  etnalyzing  in  details  the  results, 
it  seems  that,  when  aoo/c^o  increases,  the  correction  factor  reaches  an  asymptotic  value. 
Thus,  the  following  correcting  factor  has  been  chosen 


f{aoo/a[)  =  1  +  3.9  [1  -  exp(-0.035  (aco/a^  -  1))]  (4-8) 


Figure  15  presents  the  corrected  reduced  droplet  lifetime  for  the  various  computational 
cases,  T*^rrected  =  ^o/i^o  / o^q))-  This  lifetime  may  now  be  easily  correlated  to  the 

reduced  critical  temperature  by  means  of  a  linear  function  (represented  with  the  dashed 
line  in  the  figure),  the  expression  of  the  droplet  lifetime  is  given  by 

r  =  [0.0115  +  0.542  (1  - 1^)]  ^  /(a^/aE,)  (4.9) 

“o 

The  average  error  from  this  relation  is  5.2  %,  while  the  mciximum  error  for  the  presented 
results  is  around  10  %.  Application  of  Eq.(4.9)  constitutes  then  an  excellent  estimation 
technique  of  the  droplet  lifetime  for  supercriticeJ  regimes.  In  addition,  it  may  be  extended 
to  moderate  supercritical  pressures  (filled  dots  of  Fig.  15)  with  satisfying  agreement. 
The  expression  of  T*  may  be  related  to  a  Spalding  number,  which  has  more  physical 
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Oxygen  Hydrogen  Water 

Pc  {atm)  49.8  12.8  218.3 

Tc  {K)  154.6  33.1  647.3 

Table  1,  Critical  properties 


significance,  as  follows 

Tc  — To  1  —  r* 

This  correlation  shows  clearly  that  pressure  affects  the  droplet  lifetime  through  its  effect 
on  the  mixture  critical  temperature  and  on  the  ambient  thermal  diffusivity. 

5.  Vaporization  in  Hydrogen/ Water  Mixtures 
In  a  reel!  rocket  engine  combustion  chamber,  the  ambient  environment  seen  by  a  liquid 
oxygen  droplet  is  not  exclusively  composed  of  hydrogen.  One  must  consider  at  least  the 
main  combustion  product  issued  from  the  reactive  mechanism,  that  is  water.  However 
the  resolution  of  the  problem  is  not  strmghtforward  and  needs  further  considerations. 

5.1.  Condensation  of  Water 

Water  critical  properties  are  much  higher  than  the  ones  of  both  hydrogen  and  oxygen 
(See  table  5.1).  Water  is  much  less  volatile  than  oxygen  and  hydrogen.  Consequently  its 
liquid  state  exists  for  higher  temperatures  and  pressures.  Applied  to  our  problem,  it  may 
have  the  following  consequence:  when  gaseous  water  diffuses  towards  the  oxygen  droplet 
surface,  a  non-negligible  amount  of  gaseous  water  may  reach  low  temperature  regions 
so  that  water  parti2Ll  pressure  becomes  bigger  than  its  saturation  pressure.  Gaseous 
water  may  then  experience  a  phase  change  leading  to  formation  of  ice/water  particles. 
There  is  no  way  to  elude  the  problem  by  aissuming  that  condensation  process  is  negligible 
because  the  incompatibility  between  gaseous  water  concentrations  and  local  temperatures 
will  lead  to  somehow  unsolvable  problems  of  density  computations.  In  the  following 
analysis,  we  consider  that  gaseous  water  may  experience  only  one  phase  change,  that  is 
condensation;  transition  from  liquid  to  solid  state  or  direct  transition  from  gaseous  to 
solid  state  are  not  considered. 
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The  condensation  mechanism  is  complex  and  involves  different  processes  such  as  nu- 
cleation  and  particles  growth  by  germination.  In  addition,  particles  may  interact  with 
the  gaseous  flow  and  with  other  particles:  particles  coalescence.  The  purpose  here  is  not 
to  treat  the  entire  problem  induced  by  condensation  because  it  will  require  the  develop¬ 
ment  of  a  sophisticate  spray  model  which  may  overlap  our  initial  study  set.  However, 
by  analyzing  the  features  of  all  the  processes  involved,  a  simplified  modelization  of  the 
condensation  Ccin  be  obtained. 

Particle  formation  is  attributed  to  nucleation  of  molecular  clusters  called  embryos. 
As  the  saturation  ratio,  defined  as  the  ratio  of  the  partial  water  vapor  pressure  to  the 
water  saturation  pressure  {S  =  PH2o/Psat)  increases  above  unity,  embryos  larger  than  a 
critical  size  become  stable  droplets,  see  Reist  (1984).  Droplets  of  critical  size,  Rpc  are 
called  nuclei,  and  can  be  estimated  from  the  Kelvin  equation  given  by  Me  Donald  (1962) 

=  RuPpTlaS 

O’  is  the  surface  tension  of  water  and  pp  its  density,  respectively  approximately  equal 
to  7.5. 10~^  N/m  2uid  1000  kg/m^.  The  expression  is  defined  only  for  saturation  ratios 
greater  than  one.  For  saturation  ratios  less  than  one,  embryos  disintegrate  as  rapidly  as 
they  form.  Figure  16  present  results  of  Eq.(5.1)  for  different  representative  temperatures. 
We  consider  that  nucleation  process  is  not  influence  by  foreign  molecules  (oxygen  and 
hydrogen)  and  is  thus  homogeneous;  the  nucleation  rate  per  unit  volume  is  determined 
by  statistical  thermodynamics,  see  Rogers  &  Yau  (1989) 

J  =  UI&.Z  15.2) 


V27rWJiuT  W 


ZRuT 


where  the  factor  Z  denotes  the  Zeldovich  factor,  which  has  a  numerical  value  in  order  10^ 
when  all  quMtities  are  measured  in  SI  units.  Af  is  the  Avogadro  number.  By  combining 
Eqs.(5.1)  and  (5.2),  the  nucleation  rate  is  expressed  as  a  function  of  the  saturation 
ratio  and  temperature.  We  may  express,  for  a  given  gaseous  water  density  pn^O}  the 
characteristic  time  of  nucleation  Tnucieation  defined  by  the  relation 


PH^O  —  Pp  '^nucleatio 


From  the  three  previous  relations,  this  characteristic  nucleation  time  may  be  expressed 
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as  a  function  of  temperature  and  saturation  ratio.  Figure  17  presents  for  different  tem¬ 
peratures  the  nucleation  time  function  of  saturation  ratio.  This  characteristic  time  is 
extremely  sensitive  to  the  saturation  ratio.  Thus  for  saturation  ratios  above  10  it  be¬ 
comes  smaller  than  10“^°  $  and  nucleation  process  may  be  considered  as  instantaneous 
with  respect  to  our  vaporization  process  which  characteristic  time  is  close  to  10”®  s.  For 
such  saturation  ratios,  the  critical  nuclei  radius  is  between  10”^°  m  and  10“®  m,  each 
nuclei  is  then  composed  of  a  few  water  molecules. 

When  expressing  nucleation  time  function  of  temperature  and  saturation  ratio,  we 
have  expressed  the  saturation  pressure  as  a  function  of  temperature.  Such  expres¬ 
sions  derived  from  experimental  data  are  available  in  the  literature,  see  for  example 
Reid,  Prausnitz  &  Poling  (1987).  In  these  expressions,  saturation  pressure  exhibits  an 
exponential  dependency  with  respect  to  temperature:  it  vanishes  rapidly  as  temperature 
decreases.  When  liquid  oxygen  droplet  vaporizes,  the  region  which  is  likely  to  shelter  the 
water  condensation  process  corresponds  most  of  the  time  to  a  large  temperature  gradient 
region.  Then,  considering  that  nucleation  process  occurs  for  saturation  ratios  around  10 
or  occurs  at  local  thermodynamic  equilibrium  (S  =  1)  will  generate  only  small  discrepan¬ 
cies.  In  addition,  by  only  taddng  into  account  homogeneous  condensation,  we  may  have 
overestimated  the  characteristic  nucleation  times  for  saturation  ratios  close  to  one.  As 
a  matter  of  fact  gaseous  oxygen  and  hydrogen  molecules  are  expected  to  influence  the 
nucleation  process.  Based  on  these  considerations,  we  have  modeled  very  simply  nucle¬ 
ation  process  as  occurring  for  saturation  ratios  equal  to  one,  and  almost  instantaneously. 
This  hypothesis  is  equivalent  to  consider  local  thermodynamic  equilibrium.  Some  com¬ 
putations  were  carried  out  with  nucleation  process  occurring  for  saturation  ratios  above 
1  and  demonstrated  only  negligible  differences. 

After  their  formation,  condensed  particles  interact  with  the  gaseous  flow.  Water  par¬ 
ticles  are  submitted  to  thermophoretic  and  viscous  forces.  The  former  effect  tends  to 
attract  heavier  molecules  or  small  particles  to  low  temperature  regions,  as  was  discussed 
by  Talbot  et  al  (1980).  Condensed  particles  will  thus  be  attracted  by  oxygen  droplet  sur¬ 
face.  Expressions  for  thermophoretic  forces  are  bcised  on  the  derivations  of  Epstein  (1929) 
corrected  by  Brock  (1962)  to  take  into  account  non-quiescent  environment,  the  expression 
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of  the  one-dimensional  ladial  thermophoretic  force  is  obtained 


127r/i%C,  (A  -F  Ctl^)  1  dT 
~  p{l  +  SCmPh)  (1  +  2A  +  2Ct Ah)  T  dr 


(5.4) 


where  non  subscribed  quantities  refer  to  the  gaseous  environment  and  A  represents  the 
ratio  of  gas  thermal  conductivity  to  particle  thermal  conductivity.  The  Knudsen  number 
J5z  is  defined  as  the  ratio  of  molecular  mean  free  path  to  particle  radius 


Kfi  =  —  =  — 


2/i 


(5.5) 


In  Eq.(5.4)  the  constants  C,,  Ct  and  Cm  are  derived  from  kinetic  theory  and  have  been 
assigned  the  values  of  1.17,  2.18,  and  1.14,  respectively.  Brock’s  equation  yields  good 
predictions  whatever  the  ratio  A  is.  Moreover,  Talbot  et  al.  (1980)  has  shown  that  it 
exhibits  only  3  %  error  in  the  collisionless  limit,  Ah  — ^  oo,  although  it  was  originally 
designed  for  small  values  of  the  Knudsen  number.  Since  nuclei  are  typically  of  order  of 
the  molecular  mean  free  path,  viscous  forces  were  modeled  using  the  Stokes-Cunningham 
equation 

f.  =  (.-»,)  (5.6) 


where  v  and  Vp  are  the  flow  and  particles  velocities,  respectively.  The  Cunningham 
correction  coefEcient  Cc,  defined  as 


Cc 


1  +  25i 


1.257  -b  0.4  exp 


(5.7) 


has  been  introduced  by  Carlson  &  Haglund  (1964)  into  the  Stokes-flow  expression  to 
account  for  non-continuum  effects  coincident  with  water  particles  sizes  of  the  molecular 
mean  free  path  in  the  suspending  fluid.  Assuming  that  thermophoretic  and  viscous  forces 
are  the  only  forces  eicting  on  the  particle,  the  equation  of  motion  describing  particle 
dynamics  can  be  written  as 

+  (5.8) 

Actually  a  termineJ  velocity  is  reached  when  the  net  force  acting  on  particles  is  equal  to 
zero.  The  expression  is  obtained  by  equating  Eqs.(5.4)  and  (5.6).  Let  us  now  consider 
the  time  needed  to  reach  this  mechanical  equilibrium  which  expression  may  be  eetsily 
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obtained  from  previous  equations 

-  o^ 

^mecbanical  equilibrium  “  28/i 

Figure  18  presents  this  mechanical  equilibrium  characteristic  time  as  a  function  of  particle 
radius.  Since,  once  again,  this  cheiracteristic  time  is  by  several  order  of  magnitude  smaller 
than  the  characteristic  time  of  liquid  oxygen  vaporization,  we  may  consider  that  particle 
velocity  is  almost  instantaneously  equal  to  its  terminal  velocity  given  by 

2CcCsM(A  +  CtMi)  IdT 
Vp,UTminal  V  ZCml^)  (1  +  2A  +  2CtIil)  T  Sr 

Numerical  experiments  show  that  the  second  term  on  the  right  hcind  side  of  Eq.(5.10) 

is  much  smaller  than  the  first  one.  Thus  particle  velocity  is  very  close  to  gaseous  flow 

velocity.  Furthermore,  since  the  particle  radius  effect  is  entirely  included  in  this  term,  the 

overall  effect  of  condensed  particle  size  is  then  expected  to  be  small  or  even  negligible.  A 

full  treatment  of  water  condensation  process  would  require  a  different  treatment  of  each 

condensed  particle  depending  on  its  size.  Since  each  paurticle  would  grow  or  decrease 

depending  on  the  local  saturation  conditions,  and  since  each  particle  may  collide  with 

other  particles,  tracking  each  particle  size  and  position  would  necessitate  a  sophisticate 

modelization.  From  previous  considerations  such  modelization  seems  unnecessary  and 

we  may  consider  that,  as  a  first  approximation,  each  condensed  particle  has  the  same 

radius,  taken  to  be  equal  to  10“^  /im. 

The  condensation  process  may  now  be  modeled  simply.  It  can  be  reduced  to  the 
addition  of  mass  sink  term,  —mcond,H20  in  the  balctnce  equation  of  mass,  Eq.(2.1),  and  a 
source  term,  AHyap,H20  '^cond,H20  in  the  balance  equation  of  energy,  Eq.(2.3).  The  clear 
parallelism  existing  between  water  condensation  process  and  oxygen  vaporization  process 
leads  us  to  express  these  terms  with  a  method  similar  to  the  one  described  previously  in 
Eq.(3.3) 


'>^cond,H20  =  Kcond,H20  {PH2O  “  Psai)  (S-H) 

However,  a  significant  difference  exists  between  both  processes.  Water  condensation 
process  does  not  occur  at  a  fixed  location.  Thus  Eq.(5.11)  must  be  applied  for  the  gaseous 
region  adjacent  to  the  droplet  surface  and  limited  by  the  line  above  which  saturation  ratio 
becomes  lower  to  one.  This  condensation  limit  is  determined  by  means  of  an  iterative 
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procedure.  Once  created,  condensed  particles  are  convected  by  use  of  Eq.(5.10)  and  reach 
zones  where  saturation  ratio  is  lower  than  unity.  In  this  case,  condensed  particles  are 
assumed  to  disappear  almost  instantaneously. 

Neglecting  the  contribution  of  Brownian  motion  of  condensed  particles  to  the  pressure, 
we  still  can  use  Eq.(2.8)  where  p  is  the  gas  density.  Besides,  we  assume  that  the  total 
volume  is  the  sum  of  each  phase  volume.  The  condensed  mass  volume  is  related  to 
its  mass  by  the  particle  density  which  is  considered  as  independent  of  pressure  and 
temperature.  The  numerical  implementation  of  condensation  model  requires  then  the 
insertion  of  three  other  quantities:  for  each  grid  cell,  we  compute  the  condensed  particle 
mass,  volume  and  velocity. 


5.2.  Results  and  Discussion 

We  consider  now  the  vaporization  of  liquid  oxygen  droplet  in  gaseous  mixture  of  hydrogen 
and  water  for  both  subcritical  and  supercritical  regimes  of  vaporization.  In  this  latter 
case,  since  for  ternciry  mixtures  critical  composition  line  lose  most  of  its  signification,  only 
the  critical  interface  based  on  critical  temperature  is  used.  The  condensation  process 
does  not  influence  significantly  the  vaporization  process  so  that  results  are  qualitatively 
similar  to  previous  ones.  This  phase  ch£inge  experienced  by  gaseous  water  prevents  the 
water  to  reach  oxygen  droplet  surface.  Gaseous  water  is  consumed  by  the  condensation 
process  away  from  the  droplet  surface,  and  condensed  particles  are  blown  away  by  the 
gaseous  flow.  Thus,  only  oxygen  and  hydrogen  are  present  at  the  droplet  surface  and 
the  behavior  of  the  thermodynamic  equilibrium  is  very  close  to  the  one  described  in  the 
previous  section.  In  the  Fig.  19,  instantaneous  distributions  of  compositions  are  reported 
at  different  times  for  p  =  100  atm  and  To©  =  1000  K,  The  gaseous  compositions  are 
reported  to  total  gaseous  mass  whereas  the  condensed  phase  mass  fraction  is  reported 
to  the  total  mass  of  the  system.  However,  condensed  maiss  fraction  remains  very  small 
during  the  entire  droplet  lifetime,  eventually  dropping  to  zero  as  droplet  disappears.  The 
accumulation  of  condensed  water  in  the  neighborhood  of  droplet  surface  is  also  rather 
small.  Gaseous  water  does  not  diffuse  towards  the  oxygen  droplet  surface.  For  the  same 
conditions  as  Fig.  19,  Figure  20  presents  the  instantaneous  distributions  of  temperature. 
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Despite  the  fact  that  condensation  process  liberates  energy  in  the  gaseous  phase,  no 
influence  is  noted  on  the  temperature  profiles,  showing  once  again  that  condensation 
process  is  very  weak. 

Figure  21  presents,  for  different  ambient  pressures  and  temperatures,  the  instantaneous 
location  of  the  condensation  limit.  The  radius  is  reduced  by  the  instantaneous  droplet 
radius.  Condensation  first  occurs  very  close  to  the  oxygen  droplet  surface,  for  relatively 
low  temperatures.  Then,  the  limit  moves  away  from  the  droplet  reaching  higher  temper¬ 
ature  zones.  Effect  of  ambient  pressure  on  the  location  of  the  condensation  limit  may  be 
simply  explained  as  follows.  The  condensation  region  is  characterized  by  the  equality  of 
partial  pressure  of  water  ^tnd  its  saturation  pressure.  As  the  ambient  pressure  increases, 
partial  pressure  of  water  increases  also  implying  that  the  temperature  of  the  condensation 
limit  increases.  The  condensation  limit  is  then  moved  further  from  the  oxygen  droplet 
surface.  Same  reasoning  may  be  done  to  explain  ambient  temperature  effect. 

To  conclude  with  qualitative  cispects  of  condensation  process  we  can  discuss  some 
of  the  hypotheses  of  modelization.  First  computations  were  carried  out  with  different 
particle  radii,  no  significant  changes  were  observed.  In  the  model  we  only  considered 
one  phase  change.  The  treatment  of  condensed  particle  dynamics  does  not  consider 
whether  particles  are  liquid  or  solid.  Considering  a  different  phase  change  would  affect 
the  local  thermodynamic  equilibrium  and  the  liberation  of  heat  in  the  gaseous  phase. 
Since  the  overall  aspect  of  condensation  process  has  been  demonstrated  to  be  small, 
the  additive  effect  of  ice  particle  formation  is  not  expected  to  change  significantly  the 
presented  results. 

Figure  23  presents  the  droplet  lifetime  as  a  function  of  reduced  pressure  for  three  dif¬ 
ferent  ambient  compositions,  and  three  different  ambient  temperatures.  Hydrogen  com¬ 
position  is  chosen  from  Yhtj,©©  =  0.25,  0.50,  and  0.75,  respectively.  Ambient  temperature 
is  chosen  from  Too  =  1000,  1500,  and  2000  if,  respectively.  Pressure  varies  between  10 
and  300  atm.  First,  the  curves  present  obviously  the  same  trend  as  the  oxygen/hydrogen 
case.  Second,  the  effect  of  ambient  composition  is  weak.  In  the  Fig.  24  we  have  reported 
the  reduced  droplet  lifetime  using  the  same  procedure  as  for  the  pure  hydrogen  case.  The 
lifetime  is  reduced  by  means  of  initial  liquid  thermal  diffusivity  and  initieil  droplet  radius 
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with  the  correction  factor  of  Eq.(4.8).  Since  the  straight  line  appearing  in  this  figure  in 
the  same  that  in  Fig.  15,  droplet  lifetime  may  be  estimated  by  use  of  Eq.(4.9).  Average 
and  maximum  deviations  are  respectively  4.8  %  and  12  %  for  the  presented  results. 

Since  condensation  mechanism  introduces  an  other  length  scale,  the  droplet  lifetime  is 
not  simply  proportional  to  the  initial  squared  diameter.  However  condensation  process 
has  been  proven  sufEciently  small  to  neglect  this  phenomenon. 

6.  Concluding  Remarks 

A  comprehensive  analysis  of  LOX  droplet  vaporization  in  hydrogen  and  water  environ¬ 
ments  at  both  sub-  and  supercritical  conditions  has  been  conducted.  A  special  attention 
has  been  devoted  to  transport  and  thermodynamic  properties  evaluation.  In  the  case 
where  ambient  gas  in  partly  composed  of  water,  a  condensation  model  has  been  imple¬ 
mented.  Results  of  this  model  show  that  condensation  occurs  is  the  neighborhood  of 
oxygen  droplet  surface.  The  main  implication  of  this  phenomenon  is  that  gaseous  water 
does  not  reach  oxygen  droplet  surface,  gaseous  water  forms  small  condensed  particles 
which  are  blown  away  by  the  gaseous  flow.  The  overall  effect  of  condensation  is  however 
very  smadl  and  no  agglomeration  of  condensed  particles  near  the  droplet  surface  occurs. 
Effects  of  ambient  and  initial  conditions  has  been  investigated  in  detail.  Results  indicate 
that  pressure  influence  is  different  depending  upon  whether  the  vaporization  regime  is 
sub-  or  supercritical.  In  the  case  where  supercritical  regime  happens  almost  instanta¬ 
neously,  the  droplet  lifetime  may  be  estimated  simply  by  means  of  a  correlation  involving 
ambient  and  initial  conditions.  This  correlation,  derived  from  physical  considerations, 
presents  an  average  deviation  of  about  5  %. 
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Temperature,  K 

Figure  1.  Vapor-liquid  phase  equilibrium  compositions  for  O2/H2  at  various  pressures. 


Temperature,  K 

Figure  2.  Enthalpy  of  vaporization  of  O2  in  an  equilibrium  mixture  of  O2  and  H2  at  various 

pressures. 
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Figure  3.  Time  variations  of  droplet  surface  temperature  at  various  pressures;  Too  =  1000  K; 

To  =  90  K;  Do  =  100  fim. 


Reduced  Radial  Distance,  r  /  Rq 

Figure  4.  Instantaneous  distributions  of  mixture  specific  beat  in  tbe  entire  field  at  various 
times;  To  =  90  K,  Too  =  1000  iT,  p  =  100  atm,  Do  =  100  fim. 
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To  =  90  K,  Too  =  1000  K,p  =  100  atm,  Do  =  100  pm. 
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Figure  7.  Time  variations  of  the  square  of  the  reduced  diameter  at  various  pressures;  critical 
interface  based  on  temperature;  Too  =  1000  K;  To  =  90  K;  Do  =  100  fim. 


Time,  ms 

Figure  8.  Time  variations  of  the  square  of  the  reduced  diameter  at  various  pressures;  critical 
interface  based  on  composition;  Too  =  1000  K;  To  =  90  K;  Do  =  100  fim. 
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Reduced  Pressure,  p  /  Pp 

Figure  9.  Droplet  lifetime  versus  reduced  pressure  for  various  ambient  temperatures; 

To  =  90  K,  Dl  =  100  pm. 


Reduced  Pressure,  p/  Pp 

Figure  10.  Droplet  lifetime  versus  reduced  pressure  for  various  ambient  temperatures; 
To  =  90  if;  critical  interface  based  on  composition. 
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Reduced  Pressure,  p/p<j  . 

Figure  11.  Influence  of  cross-diffusion  effects  on  =  lOM^ 

pressures;  critical  interface  based  on  temperature;  =  1000  K,To- 90  A,  Uo 


Ambient  Temperature, 


Supercritical  vaporization  of  liquid  oxygen  droplets 


Reduced  Pressure,  p  /  p- 

Figure  13.  Limit  temperature  of  occurrence  of  supercritical  regime  of^vaporization  for 
initial  droplet  temperatures;  To  =  90,  and  110  K. 
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Figure  14,  Dimensionless  droplet  lifetime  function  of  (Too  —  Tc)/(Too  —  To). 
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Figure  15.  Dimensionless  corrected  droplet  lifetime  function  of  (Tc 
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Saturation  Ratio,  S  *  Ph2o  !  Psat 

Figure  17.  Characteristic  nudeatioiL  time  function  of  saturation  ratio  for  different 

temperatures. 


Particle  Radius,  m 

Figure  18.  Characteristic  mechanical  equilibrium  time  function  of  particle  radius  for  different 

temperatures. 
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Reduced  Radial  Distance,  r/  Rq 

Figure  19.  Instantaneous  distributions  of  compositions  in  the  entire  held  for  various  times, 

LOXIH2IH2O  system. 


Reduced  Radial  Distance,  r/  Rq 

Figure  20.  Instantaneous  distributions  of  temperature  in  the  entire  field  for  various  times, 

LOXIH2IH2O  system. 
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Figure  21.  Reduced  radius  of  the  limit  condensation  position  for  different  temperatures  and 

pressures. 


Figure  22,  Temperature  of  the  limit  condensation  position  for  different  temperatures  and 

pressures. 
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Reduced  Pressure,  p/pg 

Figure  23.  Droplet  lifetime  function  of  reduced  pressure  for  different  ambient  temperatures 

and  compositions. 


Effects  of  Pressure  on  Hydrocarbon  Fuel  Droplet 
Vaporization  in  Air 

Patrick  Lafon  Vigor  Yang 


1  Introduction 

In  previous  publications  we  have  presented  a  numerical  model  for  treating  droplet 
vaporization  and  combustion  in  quiescent  environments  over  a  broad  range  of  ther¬ 
modynamic  conditions.  Both  subcritical  and  supercritical  regimes  of  vaporization  axe 
studied  in  detail.  The  numerical  scheme  appears  quite  efficient  and  allows  for  a  system¬ 
atic  investigation  into  various  underlying  mechanisms  involved  in  droplet  vaporization. 
Such  studies  were  carried  out  for  the  oxygen/hydrogen  system.  The  goal  of  the  present 
study  is  to  carry  out  the  same  kind  of  systematic  studies  for  single-component  hydro¬ 
carbon  droplets  vaporizing  in  air. 

Three  hydrocarbon  fuel  droplets  are  considered,  n-pentane,  n-heptane  and  n-decane 
droplets,  respectively.  These  three  hydrocarbon  species  represent  a  good  simple  in 
terms  of  molecular  weight.  Let  us  summarize  the  characteristics  of  interest  of  all  the 
species  involved: 

2  Results  and  Discussion 

A  series  of  computations  have  been  carried  out  for  the  three  hydrocarbons.  Different 
initial  droplet  temperatures,  different  ambient  temperatures  and  pressures  have  been 


n-pentane 

n-heptane 

n-decane 

Nitrogen 

Oxygen 

w 

72.15 

100.2 

142.3 

28.01 

32.00 

Tc  {K) 

469.8 

540.3 

617.5 

126.2 

154.6 

Pc  {atm) 

33.31 

26.94 

20.70 

33.46 

49.77 

Table  1:  Species  characteristics 
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considered.  In  the  present  study  of  pure  vaporization,  there  is  only  one  length  scale 
involved  (i.e.,  droplet  diameter).^  Thus,  a  dimensional  analysis  shows  that  the  droplet 
lifetime  is  proportional  to  the  square  of  the  initial  droplet  diameter.  Therefore,  compu¬ 
tations  were  carried  out  only  for  droplets  with  an  initial  diameter  of  100  //m.  This  does 
not  mean,  however,  that  the  <i^-law  is  valid  regardless  of  the  ambient  conditions.  This 
simply  means  that  the  transient  diffusion  processes  and  time  evolution  of  the  surface 
conditions  may  be  scaled  by  the  square  of  the  initial  droplet  diameter. 

Figures  1,  2,  and  3  present  droplet  lifetimes  as  a  function  of  reduced  pressure  for 
various  ambient  temperatures,  for  n-pentane,  n-heptane,  and  n-decane  droplets,  re¬ 
spectively.  Droplet  lifetime  presents  a  similar  sensitivity  to  ambient  pressure  regard¬ 
less  of  the  hydrocarbon  species.  The  major  point  to  be  noted  is  the  difference  between 
low  ambient  temperature  cases  and  higher  ambient  temperature  cases.  For  low  ambi¬ 
ent  temperatures,  for  example  Too  =  400  K,  the  droplet  lifetime  first  increases  with 
pressure,  then  levels  off  and  eventually  decreases  a  little  bit  for  high  pressures.  This 
behavior  remains  noticeable  although  with  less  amplitude  for  ambient  temperatures  up 
to  750  iif.  For  higher  ambient  temperatures,  the  droplet  lifetime  decreases  persistently 
with  increasing  pressure.  Eventually,  when  both  ambient  temperature  and  pressure 
are  high  enough,  the  droplet  surface  may  reach  its  critical  value.  In  Figs.  1,  2,  and  3, 
the  filled  symbols  correspond  to  cases  where  supercritical  vaporization  is  reached  dur¬ 
ing  the  droplet  lifetime.  The  dashed  line  represents  the  limit  of  occurrence  of  both 
vaporization  regimes. 

Let  us  first  focus  our  attention  onto  the  thermod3mamic  equilibrium  process  in¬ 
volved  at  the  interface  tmder  subcritical  vaporization  conditions.  The  numerical  scheme 
uses  a  kinetic-like  relation  for  the  expressions  of  the  species  mass  flow  rate.^’  ®  This 
technique  allows  a  simultaneous  treatment  of  the  thermodynamic  equilibrium  existing 
at  the  interface  under  subcritical  conditions,  and  of  the  flux  conservation  at  the  inter¬ 
face.  In  addition,  it  enables  to  treat  more  complex  systems  than  binary  systems.  In  a 
previous  published  work  ^  concerning  n-pentane  droplet  combustion  in  air,  only  a  bi¬ 
nary  system  was  considered,  this  was  justified  since  due  to  the  combustion  mechanism 
oxygen  does  not  reach  the  surface.  Here,  we  have  to  handle  a  ternary  system,  which 
thermodynamic  variance  is  three,  and  which  is  thus  fully  determined  for  a  given  tem¬ 
perature,  pressure  and  one  species  composition  in  one  of  the  two  phases.  Computations 
carried  out  with  the  ternary  system  involved  (Hydrocarbon/Nitrogen/Oxygen)  show 
that  although  oxygen  mole  fraction  in  both  liquid  and  gaseous  phase  is  small,  it  pro¬ 
duces  some  important  departures  from  the  binary  system  (Hydrocarbon/Nitrogen). 
For  example,  instantaneous  conditions  of  the  thermodynamic  equilibrium  are  moni¬ 
tored  by  surface  composition  history,  and,  contrary  to  the  binary-system  case  are  no 
longer  one-to-one  function  of  the  surface  temperature. 

In  Figs.  4  and  5,  we  have  reported  for  the  n-heptane/air  system  the  surface  composi¬ 
tion  history  as  a  function  of  surface  temperature  for  two  ambient  pressures,  respectively 
p  =  5,  200  atm.  As  mentioned  previously,  compositions  are  not  one-to-one  functions 
of  temperature  as  they  were  for  binary  systems.  Besides,  for  a  given  pressure,  critical 
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coordinates  should  not  be  unique.  However,  as  temperature  increases,  both  liquid  and 
gaseous  oxygen  compositions  at  the  interface  becomes  smaller  and  general  equilibrium 
behavior  becomes  similar  to  the  one  of  a  binary  system.  In  this  case,  a  raise  of  surface 
temperature  is  linked  with  an  increase  of  n-heptane  gaseous  composition,  that  is  a 
decrease  of  nitrogen  gaseous  composition. 

Let  us  consider  figures  6  and  7  in  which  we  have  reported  respectively  final  sur¬ 
face  temperature  and  the  corresponding  enthalpy  of  vaporization  as  a  function  of  the 
ambient  pressure.  This  mixture  enthalpy  of  vaporization  is  computed  from  the  relation: 


Si  iiii 


(1) 


Two  different  ambient  temperatures  were  considered,  respectively  Too  =  400  K  and 
Too  =  lOOOiT.  For  low  ambient  temperatures,  as  pressure  increases,  the  droplet  sur¬ 
face  temperature  first  increases  then  levels  off,  which  may  be  explained  as  follows. 
As  pressure  increases,  molecular  diffusion  becomes  less  efficient  and  nitrogen  gaseous 
composition  at  the  interface  is  reduced.  From  observation  of  equilibrium  curves,  in 
Figs  4  and  5,  a  reduce  in  nitrogen  gaseous  composition  goes  with  an  increase  of  sur¬ 
face  temperature.  Of  course,  surface  temperature  cannot  become  greater  than  ambient 
temperature.  So,  as  pressure  increases,  the  temperature  gradient  at  droplet  surface 
becomes  small.  In  the  mean  time,  droplet  surface  temperature  augmentation  only 
weakly  affects  the  value  of  vaporization  enthalpy.  In  conclusion,  when  ambient  pres¬ 
sure  increases,  heat  flux  conducted  decreases  so  that  droplet  lifetime  increases.  For 
high  ambient  pressures,  the  enthalpy  of  vaporization  may  be  sufficiently  reduced  so 
that  droplet  lifetime  slightly  decreases. 

For  higher  ambient  temperatures,  as  pressure  increases,  droplet  surface  temperature 
increases.  However,  since  difference  between  surface  temperature  and  ambient  temper¬ 
ature  is  important,  heat  flux  to  droplet  surface  is  always  efficient.  Furthermore,  as 
pressure  increases,  vaporization  enthalpy  reported  in  Fig.  7  shows  a  net  decrease  even¬ 
tually  dropping  to  zero  when  supercritical  vaporization  regime  is  reached.  Then,  it 
becomes  obvious  that,  for  high  ambient  temperatures,  an  increase  in  ambient  pressure 
is  followed  by  a  decrease  in  droplet  lifetime. 


3  Correlations  of  Results 

From  observations  of  Figs.  1,  2,  and  3  and  previous  discussion,  it  seems  difficult  to 
correlate  droplet  lifetime  with  pressure  regardless  of  the  ambient  temperature.  In  the 
present  section,  we  only  consider  ambient  temperatures  above  750  K.  The  basic  idea  is 
first  to  reduce  the  droplet  lifetime  based  on  the  low  pressure  value  (p  =  1  oim),  and  then 
to  correlate  the  reduced  droplet  lifetime  by  means  of  a  function  of  ambient  pressure. 
Computations  have  been  carried  out  for  different  droplet  initial  temperatures,  which 
are  however  different  for  the  three  hydrocarbon  droplets  due  to  the  different  normal 
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boiling  temperatures  of  the  three  hydrocarbon  species.  Thus,  for  n-pentane  droplets, 
To  =  275,  300  K  have  been  considered.  For  n-heptane  and  n-decane  droplets,  To  =  300, 
325,  and  350  K  have  been  considered. 

Let  us  first  compare  our  numerical  results  obtained  for  low  pressures  with  the 
predictions  of  the  quasi-steady  theory.  This  theory  gives  two  expressions  for  droplet 
lifetime,  respectively: 

T  =  £  ^[ln(l  +  Bt)T'  (2) 

and 

(3) 

The  two  Spalding  numbers  involved  are  expressed  by: 


Bt  = 


Cp(roo  —  Ts) 

AH 


and 


By 


Yf, 

I-Ff. 


(4) 


By  equating  Eqs.(  2)  and  (3),  and  using  Raoult’s  law  for  linking  interface  equilibrium 
composition  and  temperature,  an  estimation  of  droplet  lifetime  may  be  obtained.  The 
quaisi-steady  theory  requires  that  the  thermal  diffusivity,  molecular  diffusion  coefficient 
and  specific  heat  are  time  and  space-constant.  These  properties  axe  estimated  by 
means  of  an  average  between  surface  and  ambient  conditions.  In  Fig.  8,  we  have 
reported  the  ratio  of  numerical  prediction  to  the  corresponding  theoretical  value.  The 
agreement  is  only  weak,  it  only  gives  the  order  of  magnitude.  This  disagreement  may 
be  explained  by  the  simplifying  hypotheses  used  in  the  quasi-steady  theory:  constant 
liquid  temperature,  space  constant  transport  and  thermodynamic  properties,  idealized 
thermodynamic  equilibrium.  As  can  be  seen  from  Fig.  8,  the  departure  from  the  quasi¬ 
steady  theory  is  maximum  for  low  ambient  temperature,  where  the  droplet  heat-up 
process  is  expected  to  occupy  a  longer  portion  of  the  droplet  lifetime.  This  heat-up 
process  is  not  taken  into  account  in  the  quasi-steady  theory. 

In  Tables  2,  3,  and  4,  droplet  lifetimes  at  low  pressure  conditions  (p  =  1  atm)  issued 
from  our  numerical  simulations  are  reported.  These  tables  may  be  used  to  extrapolate 
low  pressure  droplet  lifetimes  to  other  ambient  or  initial  conditions. 

Only  a  few  relevant  experimental  data  are  available  in  the  literature,  although 
many  experimental  works  have  been  carried  out  on  droplet  vaporization.  In  quiescent 
environment,  we  can  mention  the  work  performed  by,  for  example,  Faeth  et  al.^,  Kadota 
and  Hiroyasu4,  and  more  recently,  Hartfield^  and  Sato^.  The  work  of  Faeth  et  al  was 
only  concerned  by  the  combustion  of  binary-fuel  droplets,  whereas  the  work  of  Kadota 
and  Hiroyasu,  although  applying  to  pure  vaporization  processes,  was  conducted  under 
normal  gravity  conditions.  Considering  pure  vaporization  conditions  leads  to  droplet 
lifetime  of  the  order  of  10  to  30  seconds  for  the  droplet  size  handled  experimentally. 
On  the  other  hand,  micro-gravity  test  facilities  provide  micro-gravity  conditions  for 
a  much  shorter  amount  of  time.  Nevertheless,  Hartfield^  and  Sato^  have  conducted 
droplet  vaporization  experiments  under  micro-gravity  conditions.  Both  used  n-heptane 
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Too,  K  750  1000  1500  2000  2500 

To  =  275K  41.1  30.8  21.8  17.6  15.1 

To  =  SOOK  38.4  29.0  20.8  16.9  14.6 


Table  2:  Low  pressure  vaporization  lifetime  (ms);  n-pentane/air  system;  p  =  1  atm. 


Too,K 

750 

1000 

1500 

2000 

2500 

To  =  mK 

44.8 

32.6 

23.2 

18.7 

16.1 

To  =  Z25K 

42.2 

31.3 

22.2 

18.0 

15.5 

To  =  mK 

39.3 

29.5 

21.1 

17.2 

14.8 

Table  3:  Low  pressure  vaporization  lifetime  (ms);  n-heptane/air  system;  p  =  1  atm. 

droplets  vaporizing  in  a  nitrogen  atmosphere.  They  used  a  suspended  droplet  technique 
which  enable  them  to  stabilize  spatially  the  droplet.  Comparisons  of  our  numerical 
simulations  to  both  experimental  works  only  provide  poor  agreements.  Numerical  and 
experimental  results  only  agree  on  the  order  of  magnitude  of  the  droplet  lifetime.  This 
does  not  question  our  numerical  results  however.  First,  both  experimental  works  lead 
to  different  droplet  lifetimes  or  vaporization  constants,  i.e.  the  slope  in  the  later  stage 
of  the  droplet  lifetime  of  the  square  of  the  droplet  diameter  as  a  function  of  the  time. 
Two  sets  of  experimental  results  presented  by  Sato  show  some  disagreements  between 
each  others  of  the  order  of  a  few  tens  of  percents.  First  the  experimental  technique  used 
by  Hartfield  leads  to  an  over-shoot  of  the  pressure  and  a  non-zero  convective  velocity  in 
the  neighborhood  of  the  droplet  for  a  time  comparable  to  the  droplet  lifetime.  Second, 
for  both  experimental  studies,  some  uncertainties  remain  in  the  value  of  the  initial 
droplet  temperature,  which  is  strongly  affected  by  the  ambient  temperature  because 
of  the  droplet  formation  technique  used.  Moreover,  since  the  droplet  is  hold  by  a 
fiber,  this  fiber  tends  to  heat-up  the  interior  of  the  droplet  especially  if  the  ambient 
temperature  is  high.  This  leads  to  significant  difference  between  numerical  predictions 
and  experimental  results.  Thus,  the  best  agreements  between  numerical  predictions  and 
experimental  data  are  obtained  for  the  low-temperature  results  presented  by  Sato^.  In 
this  case,  the  discrepancy  between  numerical  predictions  and  experimental  results  in  of 
the  order  of  30  %.  It  is  however  very  important  to  mention  that  our  numerical  results 
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Too,K 

750 

1000 

1500 

2000 

2500 

To  =  300iir 

51.1 

36.3 

25.1 

20.1 

17.2 

To  =  325^^ 

48.9 

34.9 

24.2 

19.5 

16.7 

To  =  350A: 

46.5 

33.5 

23.4 

18.9 

16.2 

Table  4:  Low  pressure  vaporization  lifetime  (ms);  n-decane/air  system;  p  =  1  atm. 

present  the  same  trend  that  the  experimental  results  presented  by  Sato,  in  terms  of 
the  influence  of  bith  ambient  pressure  and  temperature  on  the  droplet  lifetime. 

Figure  9  presents  reduced  droplet  lifetime,  Tr  =  rfjp  =  i  atm,  as  a  function  of  the 
ambient  pressure  for  each  hydrocarbon.  All  the  different  ambient  temperatures  and 
initial  droplet  temperatures  reported  in  the  previous  tables  have  been  considered.  Re¬ 
sults  collapse  extraordinary  well,  so  that  reduced  droplet  lifetime  may  be  expressed  as 
a  function  of  pressure: 

Tr  =  exp  [-a  {p  -  Prej)/Pref]  (5) 

with  pref  =  1  atm.  Values  of  a,  average  and  maximum  deviation  are  summarize  in 
Table  5. 


a 

n-pentane  0.140 
n-heptane  0.084 
n-decane  0.050 


average  deviation 

0.9% 

0.6% 

0.7% 


maximum  deviation 

11.3  % 

5.7  % 

11.6  % 


Table  5:  Correlations  of  results. 


4  Concluding  Remarks 

A  comprehensive  study  concerning  hydrocarbon  droplets  vaporizing  in  air  has  been 
carried  out.  Only  pure  ^-heptane  and  pure  n-decane  droplets  have  been  considered. 
The  effect  of  pressure  has  been  estimated  and  explained.  Two  different  behavior  ap¬ 
pear  depending  on  ambient  temperature.  In  the  case  where  ambient  temperature  is 
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above  750  droplet  lifetime  may  be  correlated  for  various  initial  droplet  temperature 
and  ambient  temperature.  This  correlation  uses  estimation  of  droplet  lifetime  at  low 
pressure  and  a  simple  function  of  ambient  pressure.  For  the  presented  computations, 
the  average  deviation  between  numerical  results  and  the  correlation  is  less  than  2  %. 
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Figure  1;  Droplet  lifetime  as  a  function  of  the  reduced  pressure  for  different  ambient 
temperatures;  n-pentane/air  system,  To  =  300  K]  do  =  100  fim. 


Figure  2:  Droplet  lifetime  as  a  function  of  the  reduced  pressure  for  different  ambient 
temperatures;  n-heptane/air  system,  To  =  300  K\  do  =  100  fim. 
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Figure  3:  Droplet  lifetime  as  a  function  of  the  reduced  pressure  for  different  ambient 
temperatures;  n-decane/air  system,  To  =  300  K;  do  =  100  fim. 
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Figure  5:  Equilibrium  compositions  as  a  function  of  temperature;  n-heptane/air  sys¬ 
tem;  Too  =  1000  Ki  To  =  300  K,  p  =  200  atm. 


Pressure,  atm 

Figure  6:  Final  surface  temperature  function  of  ambient  pressure;  n-heptane/air  sys¬ 
tem;  To  =  300  . 
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Pressure,  atm 

Figure  7;  Final  mixture  vaporization  enthalpy  function  of  ambient  pressure;  n- 
heptane/air  system;  To  =  300  . 


Ambient  Temperature,  K 

Figure  8:  Comparison  of  numerical  results  with  theoretical  predictions  for  low  ambient 
pressure;  To  =  300  . 
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Reduced  Droplet  Lifetime,  t  /  t 


Figure  9:  Reduced  droplet  lifetime  function  of  ambient  pressure. 


12 


Droplet  Vaporization  at  Thermodynamic  Non- Equilibrium 

Conditions 


Patrick  Lafon  and  Vigor  Yang 
The  Pennsylvania  State  University 
Depsirtment  of  Mechanical  Engineering 

Abstract 

The  purpose  of  this  study  is  to  estimate  the  influences  of  thermodynamic  non¬ 
equilibrium  conditions  on  droplet  vaporization  processes.  The  study  focuses  on 
the  oxygen/hydrogen  system,  and  is  conducted  in  two  steps.  First,  a  droplet  va¬ 
porization  model  based  on  the  thermodynamic  equilibrium  assumption  is  used; 
its  predictions  of  the  vaporized  mass  flow  rate  are  compared  with  the  molec¬ 
ular  fluxes  issued  from  molecular  dynamics  theory.  Except  for  small  droplet 
diameters,  molecular  dynamics  does  not  appear  to  be  a  limiting  factor  to  the 
achievement  of  thermodynamic  equilibrium.  However,  since  some  uncertainties 
remains  in  the  expression  of  the  molecular  fluxes,  it  may  be  interested  to  com¬ 
pute  the  features  of  the  superheated  liquid  phase  resulting  from  non-equilibrium 
conditions.  This  is  the  purpose  of  the  second  part  of  this  study.  Thus,  the  limit 
of  existence  of  the  metastable  liquid  phase  are  estimated,  i.e.  the  thermodynamic 
stability  limit  (spinodal  limit)  and  the  kinetic  limit.  The  latter  limit  involves  the 
determination  of  the  homogeneous  nucleation  rate.  For  both  limit  conditions  of 
superheated  liquid,  the  change  in  terms  of  enthalpy  of  vaporization  with  respect 
to  the  thermodynamic-equilibrium  model  are  estimated  giving  some  insights  on 
the  influence  of  the  non-equilibrium  conditions  on  the  droplet  vaporization  pro- 
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cesses. 


1  Introduction 

During  the  last  decades,  many  work  were  devoted  to  the  understanding  of  the  droplet 
vaporization  process  under  high  pressure  conditions.  Highly  sophisticated  numerical 
models  were  developed  in  order  to  simulate  the  process.  The  need  of  numerical  models 
was  due  to  the  unsteady  features  of  the  overall  process,  to  the  importance  of  the  space 
and  time  variations  of  the  transport  properties,  as  long  as  to  the  importance  of  the 
dissolution  phenomenon  as  pressure  increases.  Taking  into  account  this  dissolution 
phenomenon  necessitates  the  use  of  a  real-gas  thermodynamic  modelization.  However, 
in  all  the  former  studies,  the  hypothesis  of  thermodynamic  equilibrium  at  the  droplet 
interface  is  assumed.  Thus,  it  is  postulated  that  the  molecular  transfers  are  fast  enough 
to  ensure  that  thermodynamic  equilibrium  is  reached.  In  other  words,  molecular  dy¬ 
namic  is  not  a  limiting  factor  in  the  achievement  of  thermodynamic  equilibrium.  The 
justification  lies  implicitly  in  the  fact  that  the  characteristic  time  of  completion  of 
the  equilibrium  is  short,  at  least  by  several  orders  of  magnitude,  compared  with  the 
characteristic  diffusion  times  ruling  droplet  vaporization. 

The  aim  of  this  work  is  to  possibly  provide  a  justification  to  this  assumption,  or 
to  estimate  the  consequences  of  any  departure  from  the  equilibrium  conditions.  The 
work  is  focused  on  the  system  oxygen  hydrogen,  but  the  same  developments  may  be 
carried  out  for  hydrocarbon  air  systems  and,  at  least,  some  of  the  conclusions  may 
be  used.  Therefore,  we  proceed  in  two  steps.  First,  we  determine  whether  or  not  the 
thermodynamical  equilibrium  can  be  reached  on  a  molecular  viewpoint,  by  looking 
closely  at  the  moleculax  transfers  at  the  interface.  Since  some  of  the  parameters  used 
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in  the  estimation  of  the  molecular  fluxes  are  only  roughly  known,  as  we  will  discuss 
later  on,  this  first  analysis  does  not  enable  us  to  conclude  definitely  on  this  problem. 
It  appears  thus  necessary  to  conduct  a  further  study. 

The  aim  of  this  second  study  is  to  estimate  what  the  limits  of  non-equilibrium 
conditions  are  both  on  a  thermodynamic  stability  and  kinetic  outlooks.  Here,  our  only 
concern  is  the  limit  of  superheat  for  the  liquid  phase.  The  thermodynamic  stability 
limit  or  spinodal  limit  involves  the  stability  of  the  liquid  metastable  phase,  and  is 
fully  determined  by  means  of  the  thermodynamic  modelization.  The  kinetic  limit  of 
superheat  involves  both  molecular  dynamics  theory  and  thermodynamic  modelization 
through  computation  of  the  nucleation  rate.  Although  the  molecular  fluxes  are  subject 
to  the  same  limitations  as  previously,  their  influence  on  the  estimation  of  the  nucleation 
rate  is  small  compared  to  the  one  of  the  thermodynamic  modelization,  so  that  a  good 
confidence  can  be  obtained  for  the  issued  predictions. 

The  changes  in  the  conditions  seen  by  the  interface  are  then  estimated  for  these 
two  limit  conditions,  spinodal  and  kinetic  limit,  respectively.  The  influence  of  non¬ 
equilibrium  conditions  of  the  vaporized  mass  flow  rate  and  thus  on  the  droplet  lifetime 
is  determined. 

2  Molecular  Fluxes  at  the  Interface 

If  the  interface  between  two  phases  is  studied  at  the  molecular  level,  it  appears  no  longer 
as  an  abrupt  discontinuity,  but  it  is  rather  the  locus  of  a  smooth  transition  between  the 
two  phases.  This  feature  is  presented  in  Fig.  1.  Since  this  transition  occurs  within  a 
few  molecular  mean  free  paths,  or  intermolecular  distance,  on  a  macroscopic  viewpoint 
the  interface  is  still  best  modeled  by  a  discontinuity  between  the  two  bulk  phases.  Our 
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concern  here  is  to  estimate  the  flux  of  molecules  crossing  the  interface  going  from  one 
bulk  phase  to  an  other,  in  one  or  the  other  way,  i.e.  condensation  or  vaporization 
flux.  The  net  flux  of  vaporization  is  then  the  difference  between  the  two  fluxes.  If 
Boltzmann’s  statistics  applies,  it  can  be  shown  that  the  molecular  rate  of  species  i  per 
unit  area  for  a  fixed-location  interface  may  be  expressed  by^ 


where  T ,  ks,  Ni,  mi  are  the  temperature,  the  Boltzmann’s  constant,  the  number  density 
of  species  i,  and  the  mass  of  molecule  i,  respectively.  After  rearrangements,  this  leads 
to  an  equivalent  expression  for  y'jv; 


jNi  = 


pXi/E^Ty^^ 
W  \2TtWi) 


(2) 


where  p,  Ru,  Xi  are  the  Avogadro  number,  the  density  of  the  mixture,  the  universal 
gas  constant,  and  the  molar  fraction  of  species  i  in  the  mixture,  respectively.  This 
expression  does  not  involve  in  its  current  formulation  any  assumption  concerning  the 
behavior  of  the  gas,  i.e.  the  ideal-gas  assumption  is  not  used.  Nevertheless,  when 
considering  a  droplet  vaporization  model,  one  must  consider  the  interface  regression 
velocity.  Schrage^®  has  shown  that  the  previous  expressions  have  to  be  corrected  by  a 
factor  r  function  of  the  ratio  Vaf  ^2RuT jW .  This  correction  factor  is  significant  only 
when  Valyj2RuTIW  >  0.01.  For  the  oxygen /hydrogen  system,  the  correction  factor 
will  only  be  significant  for  surface  regression  velocities  over  10  m/s.  Thus,  it  can  be 
neglected  in  a  first  approximation. 

Concerning  the  validity  of  expression  1,  two  major  restrictions  have  to  be  made. 
First,  the  expression  of  is  obtained  from  integration  of  the  velocity  distribution 
function  under  the  assumption  that  the  molecules  do  not  experience  any  collision.  For 
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a  dense  fluid,  this  assumption  is  certainly  not  valid,  but  one  could  argue  that  the 
number  of  molecules  initially  travelling  on  one  way  and  which  will  be  reflected  in  the 
other  direction  is  roughly  balanced  by  the  number  of  molecules  initially  travelling  on  the 
opposite  way  and  which  experience  also  a  collision.  Second,  for  similar  reasons,  if  the 
expression  of  j^i  is  to  be  used  to  estimate  the  phase-change  rate,  it  has  to  be  corrected 
by  a  so-called  accommodation  coefficient  which  takes  into  account  the  fact  that  only 
a  portion  of  the  molecules  striking  the  interface  is  actually  due  to  condensation  or 
vaporization.  The  remaining  fraction  is  due  to  reflection  of  molecules.  The  value  of 
the  accommodation  coefficient  depends  on  the  local  interface  conditions  and  therefore 
depends  on  the  system  involved.  Values  of  the  accommodation  coefficient  between 
0.5  and  1  have  been  reported  in  the  literature.^  When  experiments  show  that  the 
accommodation  coefficient  exhibits  values  very  different  from  unity,  it  is  most  likely 
that  the  molecular  dynamics  theory  used  does  not  apply  for  the  system  considered. 

A  way  to  assess  the  validity  of  the  thermodynamic-equilibrium  assumption  is  to 
compare  the  respective  orders  of  magnitude  of  the  net  vaporization  mass  flow  rate  under 
equilibrium  conditions,  and  the  mass  flow  rates,  vaporization  and  condensation,  issued 
from  molecular  dynamics  theory.  To  overcome  the  difficulty  related  to  the  uncertainty 
associated  with  the  accommodation  coefficient,  we  only  consider  the  following  ratios 

{rriijNi  {rriijNi 

where  As,  rrn,  thi^eg  represent  the  area  of  the  interface,  the  mass  of  molecule  i,  and  the 
equilibrium  net  vaporization  mass  flow  rate,  respectively.  (rriijNi  As)‘  and  {mijN-  AsY 
represent  the  mass  flow  rate  across  the  interface  coming  from  the  bulk  liquid  phase,  i.e. 
vaporization  mass  flow  rate,  and  the  mass  flow  rate  across  the  interface  coming  from  the 
bulk  gaseous  phase,  respectively.  The  estimations  of  the  vaporization  mass  flow  rate 
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under  thermodynamic-equilibrium  conditions  are  obtained  from  numerical  simulations 
issued  from  a  one-dimensional  model.  The  model  and  the  theoretical  formulation  used 
have  been  presented  by  the  authors  in  a  previous  paper.®  The  model  ha.s  been  applied 
successfully  to  the  simulation  of  oxygen  droplet  vaporization  in  hydrogen  atmospheres 
under  a  broad  range  of  ambient  pressure  conditions. 

From  Eq.  (2),  we  see  that  molecular  fluxes  are  fully  determined  by  the  bulk  phase 
temperature,  density  and  composition.  In  Fig.  2,  the  decimal  logarithm  of  the  ratios 
of  Eq.  (3)  axe  reported  throughout  the  droplet  lifetime  for  the  two  species  involved 
in  this  particular  system.  The  ambient  pressure,  the  ambient  temperature,  and  the 
initial  droplet  diameter  are  equal  to  30  atm,  1000  K,  and  100  /xm,  respectively.  The 
divergence  exhibited  by  the  two  curves  concerning  hydrogen  rates  comes  from  the  sign 
change  of  hydrogen  mciss  flow  rate.  The  hydrogen  starts  to  condensate  and  aggre¬ 
gates  within  the  liquid  droplet,  after  a  while,  the  trends  is  reversed  and  the  hydrogen 
vaporizes.  However,  for  hydrogen,  the  predictions  of  mass  flow  rate  from  molecular 
dynamics  theory  are  by  several  orders  of  magnitude  bigger  than  the  hydrogen  equilib¬ 
rium  mass  flow  rate.  Despite  the  uncertainty  in  the  accommodation  coefficient,  this 
demonstrates  that  hydrogen  transfer  at  the  interface  is  by  no  means  a  limiting  factor 
in  the  achievement  of  the  thermodynamic  equilibrium. 

For  the  case  reported  in  Fig.  2,  the  molecular  dynamics  gives  estimations  of  the  oxy¬ 
gen  mass  flow  rate,  condensation  or  vaporization,  by  slightly  more  than  3  orders  of  mag¬ 
nitude  bigger  than  the  net  vaporization  mass  flow  rate  obtained  under  thermodynamic- 
equilibrium  conditions.  The  smaller  of  the  two  ratios  is  the  one  concerning  oxygen  flux 
from  the  gaseous  phase  towards  the  liquid  phase,  i.e.  the  condensation  mass  flow  rate. 
This  features  is  essentially  due  to  the  change  in  density  between  liquid  and  gaseous 
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phase. 


For  subcritical  regimes  of  vaporization,  i.e.  when  thermodynamic  equilibrium  exists 
throughout  the  droplet  lifetime,  the  droplet  vaporization  process  is  almost  quasi-steady 
except  at  the  beginning  of  the  process.  In  this  case,  the  famous  cP.iaw  holds,  the 
droplet  surface  regression  is  linear  with  time  and  it  can  be  shown  that  the  mass  flow 
rate  is  proportional  to  the  droplet  diameter.  If  this  quasi-steady  regime,  the  surface 
temperature  reaches  a  constant  value  after  the  heat-up  period,  so  does  the  interface 
composition.  Thus,  the  molecular  dynamics  mass  flow  rate  is  proportional  to  the 
instantaneous  droplet  diameter.  Furthermore,  in  this  pure  vaporization  case,  there  is 
only  one  length  scale  involved  in  the  overall  process.^  A  dimensional  analysis  shows 
that  all  characteristic  times  are  related  to  the  squared  initial  droplet  diameter.  Thus, 
the  droplet  behavior  throughout  its  lifetime  may  be  related  to  its  initial  diameter.  From 
the  previous  remarks,  we  see  that  the  order  of  magnitude  of  the  ratios  defined  in  Eq.  (3) 
will  be  reduced  by  1  as  the  initial  droplet  diameter  is  divided  by  10.  Thus,  for  the  same 
temperature  and  pressure  conditions  as  the  ones  reported  in  Fig.  2,  the  difference  in 
order  of  magnitude  will  drop  to  2  for  an  initial  droplet  diameter  of  10  nm.  Thus,  for 
small  droplet,  the  thermodynamic  equilibrium  assumption  becomes  questionable. 

In  Figs.  3  and  4,  the  influence  of  ambient  pressure  and  ambient  temperature  is 
analyzed.  In  those  plots,  only  the  minimum  of  the  two  ratios  defined  in  Eq.  (3)  for 
oxygen  are  reported.  The  ambient  pressure  ranges  from  10  atm  to  50  atm,  and  the 
ambient  temperature  from  1000  K  to  2500  K.  The  order  of  magnitude  of  the  ratio 
between  the  two  mass  flow  rates  is  only  slightly  affected  by  the  change  in  ambient 
conditions.  The  estimations  of  the  mass  flow  rate  issued  from  the  molecular  dynamics 
theory  are  sensitive  to  the  ambient  pressure  by  the  change  induced  in  the  density,  in  ad- 
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dition,  when  the  pressure  varies,  the  equilibrium  compositions  and  temperature  of  the 
interface  are  affected.  On  the  other  hand,  a  change  in  pressure  does  not  affect  consid¬ 
erably  the  value  of  the  equilibrium  net  vaporization  rate.  Since,  as  pressure  increases, 
the  partial  density  of  oxygen  on  the  gaseous  phase  of  the  interface  increases  also,  the 
thermodynamic  equilibrium  becomes  more  and  more  valid  as  pressure  increases.  The 
ambient  temperature  as  only  a  small  influence  upon  the  interface  conditions.  However, 
its  influence  is  more  significant  on  the  droplet  lifetime  and  therefore  on  the  equilibrium 
net  vaporization  mass  flow  rate,  which  increases  significantly  as  ambient  temperature 
increases.  This  explains  why  the  order  of  magnitude  of  the  ratio  of  the  two  mass  flow 
rates  drops  slightly  with  increasing  ambient  temperatures. 

As  a  conclusion  to  this  first  part,  after  estimating  the  eflBciency  of  the  molecular 
transfer  at  the  interface,  the  thermodynamic  equilibrium  assumption  is  valid  for  oxygen 
droplets  vaporizing  in  hydrogen  atmospheres  as  long  as  the  droplet  initial  diameter  is 
above  a  threshold  value  close  to  10  fim.  Usually  the  oxygen  droplet  diameter  atomized 
in  rocket  engines  combustion  chamber  are  close  to  100  nm.  Small  droplets,  such  as 
the  ones  mentioned  previously,  are  obtained  when  99.9  %  of  the  droplet  mass  have 
vaporized,  and  such  small  droplets  are  usually  neglected.  The  achievement  of  thermo¬ 
dynamic  equilibrium  depends  on  the  ambient  conditions  of  temperature  and  pressure, 
also  this  sensitivity  is  less  significant  than  the  sensitivity  to  the  droplet  diameter.  In  this 
first  part,  we  have  considered  that  thermodynamic  equilibrium  assumption  for  droplet 
vaporization  modelization  holds  as  long  as  the  molecular  transfer  is  by  2  or  3  orders 
of  magnitude  bigger  that  the  actual  transfer  required  under  equilibrium  conditions. 
Despite  our  efforts,  some  uncertainties  still  remain  for  the  estimation  of  the  molecular 
dynamics  fluxes.  First,  there  is  the  uncertainty  in  the  validity  in  the  chosen  statistics, 


8 


i.e.  Boltzmann’s  statistics  versus  Fermi  or  Dirac’s  statistics.  Second,  there  is  the  un¬ 
certainty  in  the  order  of  magnitude  of  the  accommodation  coefficient.  It  is  believed 
although  not  proven  that  this  will  not  alFect  the  previous  conclusions.  Nevertheless,  it 
might  be  interesting  to  estimate  what  will  be  the  effects  on  the  surface  conditions,  in 
terms  of  compositions  and  enthalpy  of  vaporization,  under  non-equilibrium  conditions. 
This  will  be  the  aim  of  the  next  section. 


3  Limits  of  Superheat 


Let  us  discuss  an  interesting  and  necessary  matter  about  thermodynamic  phase  transi¬ 
tions  and  phase  stability.  In  most  of  physical  cases,  phase  changes  do  not  occur  under 
thermodynamic  equilibrium  conditions^  Whether  the  liquid  is  superheated  or  the  vapor 
is  supercooled.  For  .single-component  mixture,  the  liquid  phase  is  superheated  when, 
for  a  given  pressure,  its  temperature  is  above  saturation  value,  and  vice  versa  for 
supercooled  vapor.  Such  phases  are  called  metastable  phases.  The  state  of  a  system, 
with  given  values  of  the  total  entropy  S,  of  the  number  of  molecules  N,  and  of  the 
volume  V,  satisfies  the  relation 

8U  =  ^  (4) 


where  U  represents  the  total  internal  energy  of  the  system.  This  state  is  stable  if  and 
only  if  it  is  located  on  a  local  minimum  of  the  phase  space.  From  basic  thermodynamic 
relations,  it  can  be  shown  that  a  stable  thermodynamic  state  satisfies  several  stability 
criteria.^  For  a  single-component  mixture,  the  first  criterion  is  called  the  criterion  of 
thermal  stability  and  requires  that 

Cv  >  0  (5) 
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This  criterion  is  virtually  never  violated,  and  the  second  stability  criterion  for  a  single¬ 
component  mixture,  called  the  criterion  of  mechanical  stability  becomes  a  necessary 
and  sufficient  condition,  it  requires  that 


>0 


(6) 


For  a  binary  mixture,  like  the  one  we  are  interested  in,  an  extra  criterion  of  stability 
must  be  satisfied,  it  is  called  the  criterion  of  material  stability  or  criterion  of  diifusional 
stability.  The  former  terminology  is  preferable  since  no  diffusion  processes  are  involved. 
It  is  expressed  as  follows 


(7) 


This  relation  may  be  transformed  by  using  derivatives  of  the  chemical  potential  with 
respect  to  the  number  of  moles  of  both  species 


where  n  represents  the  total  number  of  moles  of  the  system.  These  latter  derivatives 
possess  interesting  properties  of  symmetry  with  respect  to  species  1  and  2,  due  to  the 
second  Gibbs-Duhem  equation  applied  to  the  Gibbs  free  energy.  It  is  interesting  to 
tell  that  the  criterion  of  Eq.  (7)  is  satisfied  or  violated  simultaneously  for  both  species. 
As  for  a  single-component  mixture  the  criterion  of  mechanical  stability  was  violated 
before  the  criterion  of  thermal  stability,  for  a  binary  mixture,  the  criterion  of  material 
stability  is  the  first  to  be  violated. 

When  the  system  crosses  this  limit  of  thermodynamic  stability,  it  experiences  a 
sudden  phase  transition,  i.e.  the  superheated  liquid  becomes  suddenly  vapor.  This 
thermodynamic  limit  of  superheat  is  called  the  spinodal  limit.  On  a  molecular  level, 
the  sudden  phase  change  is  due  to  internal  fluctuations  within  the  liquid  which  are  no 


10 


longer  compatible  with  the  existence  of  a  liquid  phase.  Such  fluctuations  are  called 
heterophcLse  fluctuations. 


There  is  a  second  limit  to  the  existence  of  metastable  phases  which  is  called  the 
kinetic  limit,  since  it  involves  molecular  dynamics.  The  heterophase  fluctuations,  which 
have  been  reported  before,  occurs  not  only  at  the  spinodal  limit  but  also  all  over 
the  stable  phase  and  metastable  phase  domain.  Those  heterophase  fluctuations  are 
characterized  by  some  local  molecular  density  fluctuations.  Let  us  concentrate  our 
attention  on  the  case  of  a  superheated  liquid,  which  is  oh  special  interest  for  us.  For 
this  case,  heterophase  fluctuations  will  tend  to  create  bubbles  within  the  metastable 
liquid  phase.  By  computing  the  Gibbs  free  energy  of  the  nucleus  bubble,  it  can  be 
shown  that  these  nuclei  are  either  stable  or  unstable  depending  on  their  size  and  thus 
on  their  number  of  molecules.  There  exists  a  critical  nucleus  which  is  located  on  a  local 
extrem\un  of  the  Gibbs  energy  function.  When  a  critical  nucleus  gains  one]*m^ecule, 
it  will  tend  to  grown  indefinitely,  whereas  when  it  loses  one  molecule,  it  will  tend 
to  decay.  The  homogeneous  nucleation  rate  is  the  rate  of  creation  of  nucleus  of  the 
critical  size.  At  the  saturation  conditions,  the  homogeneous  nucleation  rate  is  zero, 
essentially  be<^^es  the  radius  of  the  critical  nucleus  is  infinite.  But,  as  the  level  of 


superheat  is  increased,  this  homogeneous  nucleation  rate  will  switch  progressively  from 
extremely  low  values  to  extremely  high  value.  Eventually,  its  value  will  be  high  enough 
to  physically  observe  au  instantaneous  phase  transition.  Usually,  a  threshold  value  of 
10^^  cm“®5~^  is  arbitrary  chosen  to  constitute  the  kinetic  limit  of  superheat.  Let  us  now 
introduce  the  theoretical  basis  associated  with  the  computation  of  the  homogeneous 
nucleation  rate  for  binary  systems. 

The  index  1  and  2  denotes  quantities  with  respect  to  oxygen,  and  hydrogen,  respec- 
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tively.  The  homogeneous  nucleation  rate  is  expressed  in  the  following  form 

J  =  K  ■  ^9^ 

where  K  is  a  kinetic  factor,  and  AG*  is  the  Gibbs  free  energy  of  the  critical  nucleus.  The 
Gibbs  free  energy  of  a  nucleus  AG(Ar^,  N^)  is  a  function  of  the  number  of  molecules 
of  both  species  within  the  nucleus  bubble.  The  AG(iVj,  N^)  plane  exhibits  a  saddle 
point,  which  position  corresponds  to  the  position  of  the  critical  nucleus.  It  corresponds 
to  the  lowest  passage  over  an  energy  barrier  for  the  nucleus  to  reach  a  size  above 
which  it  will  grow  to  a  macroscopirjfopl^  The  central  assumption  of  the  classical 
nucleation  theory  is  the  Gibbs’  capillarity  assumption  which  states  that  the  bubble 
can  be  described  by  using  classical  thermodynamic.  The  expression  of  the  Gibbs  free 
energy  is  given  by  a  sum  of  a  bulk  plus  a  surface  contribution 

AG(Ar*,  Nl)  =  n‘V!*  +  E  -  E  "W  +  (10) 

x=l,2  1=1,2 

where  n^,  n,-^,  and  nf,  represent  the  number  of  moles  of  species  i  in  the  total  bubble, 
in  the  bulk  volume  of  the  bubble,  and  in  the  surface  of  the  bubble,  respectively.  These 
mole  numbers  may  be  related  to  the  molecule  numbers  by  means  of  the  Avogadro  num¬ 
ber.  a  represents  the  surface  tension  of  the  bubble.  In  this  expression,  a  distinction  has 
been  made  between  the  molecules  of  the  bubble  on  the  surface  and  in  the  bulk  volume 
of  the  bubble.  However,  the  surface  of  the  bubble  is  assumed  to  be  in  equilibrium  with 
the  bulk  volume 

lif  =  ixf  (11) 

so  that  the  Eq.  (10)  may  be  simplified  to 

^G(Nl  N‘)  =  E  (>*5  -  ><!)  +  (12) 

»=1,2 
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Furthermore,  by  using  the  Gibbs-Duhem  equation  for  the  bulk  bubble  molecules 

(13) 

j=l,2 

and  the  for  the  surface  bubble  molecules 

X;  nf  dfif  +  iirR^^'da-  =  0  (14) 

.=1,2 

it  can  be  shown  that  the  system  characterizing  the  location  of  the  saddle  point  of  the 
nucleus  Gibbs  free  energy 


\  /7v|,r,p 

(15) 

(dAG\ 

(16) 

may  be  simplified  to  give  the  so-called  Kelvin  equations 

+  ^  =  0 

(17) 

+  ^  =  0 

(18) 

where  u,-  represents  the  partial  volume  of  species  i,  defined  as 

(19) 

which  satisfies  the  following  relation 

(20) 

One  important  remark  must  be  done  concerning  the  expressions  of  Eqs.  (17)  and  (18). 
The  bubble  chemical  potentials  and  /ij  involved  concern  whether  the  surface  of  the 
bubble  or  the  bulk  volume  of  the  bubble.  In  the  latter  case,  they  can  be  computed 
from  the  characteristics  of  the  bulk  volume  of  the  bubble,  i.e.  temperature,  pressure. 
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and  composition.  On  the  other  hand,  the  partial  molar  volumes  Uj  and  uj  are  related 
to  the  composition  of  the  whole  bubble,  surface  plus  bulk  volume.  They  represent  the 
increase  of  the  volume  due  to  an  increase  of  one  in  the  bubble  of  the  number  of  moles 
1  and  2,  respectively.  From  the  bulk  volume  composition  and  classical  thermodynamic 
modelization,  an  estimation  of  the  partial  molar  volume  for  the  bulk  volume  of  the 
bubble  Vi’’  may  be  obtained.  Although  only  the  molecules  of  the  bulk  volume  of  the 
bubble  contribute  to  the  volume  of  the  bubble,  vf  and  v-  are  not  equal  but  follow 


Since,  we  do  not  have  more  details  about  the  molecule  repartition  between  the  surface 
and  the  bulk  volume,  vf  cannot  be  obtained.  However,  as  the  size  of  the  nucleus 
increases,  it  can  be  shown  easily  that  the  part  of  the  molecules  forming  the  surface 
becomes  negligible  with  respect  to  the  total  number  of  molecules  within  the  bubble,  so 
that  from  relation  21,  v-  and  vf’'  become  equal.  This  assumption,  although  only  valid 
for  large  nuclei,  will  be  done  in  the  remaining  of  this  study.  From  this  assumption, 
and  use  of  Eqs.  (12),  (17),and  (18),  the  Gibbs  free  energy  of  the  critical  nucleus  may 
be  simply  expressed  as  follows 

AG*  =  (22) 

In  addition,  the  critical  nucleus  is  in  mechanical  equilibrium  with  the  surrounding 
liquid,  so  that  the  Laplace’s  equation  holds 


P 


(23) 


and  in  thermal  equilibrium  with  the  surrounding  liquid 


% 

For  a  binary  mixture,  the  saturation  ratio  is  classically  defined  as  the  ratio  of  the 
molar  fraction  of  species  2  in  the  liquid  phase  with  respect  to  the  one  under  saturation 
conditions 

c  ^2 

s  =  — —  (25) 


For  a  given  liquid  pressure,  liquid  temperature,  and  saturation  ratio,  the  features  of 
the  critical  nucleus,  i.e.  numbers  of  molecules,  radius,  pressure,  and  temperature,  are 
fully  determined  by  solving  simultaneously  Eqs.  (17),  (18),  (23),  and  (24). 

To  compute  the  Gibbs  free  energy  of  the  critical  nucleus,  given  that  the  Gibbs’ 
capillarity  applies,  one  needs  a  thermodynamic  modelization  in  order  to  calculate  the 
thermodynamic  quantities  needed.  In  the  previous  studies  concerning  droplet  vaporiza¬ 
tion  at  high  pressure,  which  we  have  carried  out,^’®  the  Soave-Redlich-Kwong  equation 
of  state  has  been  used.  Predictions  of  the  liquid  density  are  close  to  10  %  of  the  exper¬ 
imental.  However,  the  use  of  such  a  simple  equation  of  state  was  necessary  due  to  the 
need  to  calculate  anal3d;ically  first  and  second  derivatives  of  primary  thermodynamic 
variables  such  as  Gibbs  free  energy.  In  the  present  study,  we  face  the  same  constraints. 
In  addition,  it  is  better  to  be  consistent  with  previous  studies.  Thus,  the  thermody¬ 
namic  modelization  is  based  on  th(S  ^ve-Kedlich-Kwoi^  equation  of  state.  We  malce 
the  assumption  that  such  an  equation  of  state  is  able  to  predict  the  thermodynamic 
behavior  of  a  metastable  phase.  Looijmans  et  alJ  have  conducted  the  same  kind  of 
analysis  using  Soave-Redlich-Kwong  equation  of  state,  although  it  was  for  alkanes. 
When  expressing  the  pressure  in  terms  of  temperature  and  density,  it  takes  the  form^ 

pRuT  aa 


P  = 


■  + 


(26) 


{W  -bp)  W{W  +  bp) 
where  a  and  b  accounts  for  attractive  forces  and  repulsive  forces  between  molecules. 
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respectively.  Mixing  rules  to  compute  a  and  b  are  based  on  corresponding  state  principle 

N  N  N 

OCX  =  ^  ^  ^  ]  XiX.i(l  kij')yj(ii(ijOLiOLj  ,  h  =  ]  JCj bj  (27) 

»  i  i 

where  a,-,  b,  and  Q',(7’)  are  related  to  the  critical  properties  and  acentric  factor  of 
constituent  species  as  follows 

0.08664iJ„Tc, 


a;  = 


0.42748722  Tc? 
Pci 


bi  = 


Pci 


(28) 


0,(7)  =  [l  +  M  (l  -  v/tVtJ)]"  with 


M  =  0.48  +  1.574a;,-  -  0.176a;? 


All  the  thermodynamic  properties  required  in  the  analysis  come  from  derivation  of 
Eq.(26)  and  (27).  Some  details  cam  be  found  in  previous  publications.^;® 

Finally,  in  order  to  be  able  to  compute  the  Gibbs  free  energy  of  the  critical  nucleus, 
we  need  an  estimation  of  the  surface  tension  of  the  critical  bubble.  On  one  hand 
its  estimation  is  crucial  is  the  final  accuracy  of  the  Gibbs  free  energy,  the  pressure 
difference  between  the  bubble  and  its  surrounding,  and  thus  on  the  nucleation  rate, 
on  the  other  hand  not  many  data  exist  in  the  literature.  This  is  paxticulaxly  true  for 
the  oxygen/hydrogen  system.  The  surface  tension  of  a  mixture  may  be  estimated  by 
means  of  the  Macleod-Sugden  correlation,®  it  is  expressed  by 


ri/4 


=  EiPi)  [(pxdw)'  -  (pXiiwf] 

t 


(29) 


with  W  and  <t  in  g/cm^,  gfmol  and  dynfcm.  The  parameter  Pi  are  the  so-called 
parachors  of  the  components.  Even  though  is  was  originally  suggested  to  compute  then 
from  molecular  structure,  best  agreement  is  obtained  by  fitting  experimental  results. 
For  instance,  for  oxygen,  an  experimental  expression  of  the  surface  tension  gives'^ 
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with  <T°,  and  p  equal  to  3.95  dynicm,  and  1.255  respectively,  which  gives  an  evaluation 
of  the  parachor  around  30.0.  For  hydrogen,  such  experimental  correlation  just  does  not 
exist.  In  order  to  obtain  an  estimation  of  the  parachor  for  hydrogen,  we  proceed  with 
the  following  procedure.  We  estimate  the  surface  tension  of  pure  hydrogen  by  means  of 
a  theoretical  expression  based  on  the  corresponding  state  principle  et  supposedly  valid 
for  polar  liquids.®  The  expression  introduces  the  Stiel  polar  factor  X  and  presents  the 
following  shape 

^  c  =  (31) 

with(^^w)i,  Pc)  Tc^  and  u  being  the  surface  tension  in  dynfcm,  the  critical  pressure  in 
bar,  the  critical  temperature  in  K,  and  the  acentric  factor,  respectively.  Qp,  m,  and  X 
have  respectively  the  following  expressions 

Qp  =  0. 1560  +  0.365 a; -  1. 754 X- 13.57^2- 0.506  w^  +  l. 287 w  A"  (32) 
m  =  1.210  4- 0.5385a; -14.61  A- 32.07^2 -  1.656u;=^  + 22.03a;  A  (33) 
A  =  log  [Pvp(0.6r,)/Pe]  + 1.70  a; +  1.552  (34) 

This  enables  us  to  obtain  a  value  of  the  parachor  for  hydrogen  close  to  68.0  Eventually, 
we  are  now  able  to  get  an  estimation  of  the  Gibbs  free  energy  of  the  critical  nucleus 
bubble. 

Figure  6  presents  the  contour  levels  of  the  nucleus  bubble  Gibbs  free  energy  AG, 
reduced  by  ksT,  as  a.  function  of  and  N^,  i.e.^  the  number  of  molecules  in  the 
nucleus  bubble  of  oxygen  and  hydrogen,  respectively.  The  liquid  temperature,  liquid 
pressure,  and  saturation  ratio  are  equal  to  110  K,  30  atm,  and  5,  respectively.  The 
saddle  point  position  has  been  determined  by  solving  the  system  presented  previously. 
However,  for  the  plot  of  Fig.  6,  an  approximate  method  has  been  used  to  compute  the 
Gibbs  free  energy  of  the  nucleus.  For  a  given  bubble,  the  repartition  of  the  molecules 
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between  the  surface  and  the  bulk  volume  is  unknown,  partly  because  we  do  not  know 
how  to  model  the  thermodynamic  behavior  of  the  surface  molecules.  Thus,  the  bulk 
volume  composition  of  the  bubble,  and  the  expression  of  Eq.  (10),  even  for  a  given 
number  of  molecules  within  the  nucleus,  cannot  be  obtained.  Nevertheless,  the  shape 
of  the  Gibbs  free  energy  can  be  obtained  in  the  vicinity  of  the  critical  nucleus.^’ The 
bulk  volume  composition  is  assumed  to  be  equal  to  the  total  bubble  composition,  the 
surface  tension  and  the  pressure  inside  the  bubble  are  kept  constant  at  the  value  of 
the  critical  nucleus.  Thus,  though  the  nuclei  do  not  satisfy  the  Laplace’s  equation,  the 
Gibbs  free  energy  derivatives  with  respect  to  the  molecule  number  of  each  of  the  two 
species  are  equal  to  the  left-hand  side  of  the  Kelvin’s  equation  reported  in  Eqs.  (17) 
and  (18). 

In  Fig.  6,  the  saddle  point  is  located  very  close  to  the  JVj  axis,  which  indicates 
that  the  critical  nucleus  bubble  is  composed  mainly  of  hydrogen.  For  the  saturated 
conditions  considered,  the  critical  bubble  is  composed  of  17  molecules  of  oxygen  and 
110  molecules  of  hydrogen,  corresponding  to  a  nucleus  radius  of  9.5  •  10~^°  m. 

To  compute  the  nucleation  rate,  one  must  determine  the  expression  of  the  kinetic 
factor  K  in  Eq.  (9).  Developments  were  first  carried  out  by  Reiss®  and  refined  by 
Stauffer,^^  and  reported  by  Zeng.^^  An  expression  of  K  is  obtained 

K  =  MAr^ABZ  (35) 

A  is  the  surface  area  of  the  critical  nucleus  bubble,  B  is  the  average  growth  rate,  which 
is  given  by 

Q  _ _ hiJN-i _ 

jui  sin^  <l>  -f  jjvj  cos2  (j> 

and  jV,  given  by  Eq.  (1).  Z  is  the  generalized  Zel’dovich  factor  given  by  the  following 
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equations 


^  _  -^^11  cos^  +2  Di2  cos  <j>  sin  <f>  +  D22  sin^  <t> 

{Di2  —  Dll  D22)^^^ 


with 


p. ) 
■’  2\dN*dN^j) 


T,P 


(37) 


(38) 


^  is  the  angle  between  the  iVj  axis  and  the  direction  of  growth  of  the  critical  nucleus 
bubble,  its  expression  is 

'",2  ,  mV 


1/2 


where 


tan  ^  =  s  +  5^  +  4-  , 

V  JN2J 


^  _  -P22(m/m)  -  Dll 
2  Di2 


(39) 


(40) 


In  Fig.  6,  the  location  of  the  critical  nucleus  as  well  as  the  predicted  angle  ^  is 
reported,  indicating  the  direction  of  growth  of  the  critical  nucleus. 

The  nucleation  rate  appears  to  be  very  sensitive  to  the  saturation  ratio.  For  exam¬ 


ple,  its  value  will  be  changed  by  several  order  of  magnitude  as  the  saturation  ratio  is 
changed  by  a  few  percents.  Classically,  a  threshold  value  of  the  nucleation  rate  is  arbi¬ 
trarily  chosen  to  represent  the  limit  between  finite  rate  nucleation  and  instantaneous 
nucleation.  In  this  study,  the  limit  is  J  =  10^^  -g  recommended  for 

bubble  nucleation V  whereas  a  value  of  J  =  10®  is  recommended  for  droplet 

nucleation.  Despite  its  arbitrary  nature,  the  limit  value  of  the  nucleation  rate  does  not 


influence  significantly  the  final  results.  In  Fig.  5,  the  nucleation  rate  is  reported  as  a 
function  of  the  saturation  ratio  for  different  pressures,  20,  30,  40,  50  atm,  respectively. 
In  every  case,  the  plot  looks  like  a  vertical  line,  so  that  the  nucleation  rate  goes  from 
10°  to  10^^  for  a  change  of  a  few  percents  of  the  saturation  ratio.  The  threshold 


value  of  the  nucleation  rate  cannot  be  related  easily  to  a  nucleation  characteristic  time. 
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through  an  expression  like  the  following 


Enucleation  — 


Pi 


'lElnuclcus  J 


(41) 


where  pi  and  mnucUus  represent  the  density  of  the  superheat  liquid  and  the  mass  of  the 
critical  nucleus,  respectively.  The  definition  of  the  nucleation  rate  is  the  rate  of  creation 
of  critical  nuclei  per  unit  time  and  volume,  those  nuclei  are  unstable  and  tend  to  grow 
indefinitely.  Thus,  the  characteristic  time  of  micleation  is  expected  to  be  much  shorter 
than  the  predictions  issued  from  Eq  (41)  and,  at  the  chosen  kinetic  limit  of  superheat, 
it  is  expected  to  be  much  shorter  than  any  of  the  characteristic  time  involved  in  the 
droplet  vaporization  process. 

In  the  Figs.  7-10,  the  ratio,  for  both  oxygen  and  hydrogen,  between  the  enthalpy 
of  vaporization  as  a  function  of  temperature  under  superheat  conditions  -for  both  the 
spinodal  and  the  kinetic  limit  of  superheat-  and  the  enthalpy  of  vaporization  at  the 
saturation  conditions 

Ahi/Ahi,  sat  (42) 


is  reported  for  different  pressures,  20,  30,  40  and  50  atm,  respectively.  First,  as  it  may 
have  been  expected  a  priori,  the  ratio  of  the  enthalpies  of  vaporization  becomes  bigger 
at  the  spinodal  limit.  This  limit,  as  was  mentioned  previously,  constitutes  a  limit  of 
thermodynamic  stability  and  is  never  reached  in  principle.  Let  us  focus  our  attention 
on  the  kinetic  limit  of  vaporization.  The  behavior  of  the  ratio  of  the  enthalpies  of 
vaporization  depends  strongly  of  the  species  considered.  For  oxygen,  the  ratio  remains 
very  close  to  1,  so  that  the  relative  difference  is  smaller  than  10  %.  On  the  other 
hand,  for  hydrogen,  the  ratio  becomes  very  different  from  unity  and  tends  eventually 
to  diverge  when  the  saturation  temperature  of  pure  oxygen  is  reached.  Although  this 
latter  feature  seems  to  show  that  non-equilibrium  conditions  have  significant  influence 
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on  the  droplet  vaporization  process,  this  feature  must  be  balanced  by  the  fact  that,  for 
the  oxygen/hydrogen  system,  the  ratio  is  always  negligible,  as  might  have  been 

deduced  from  Fig.  2.  In  fact,  only  the  ratio  of  enthalpies  of  vaporization  concerning 
oxygen  should  be  considered.  In  this  case,  non-equilibrium  conditions  do  not  affect 
by  more  than  10  %  the  enthalpy  of  vaporization  at  the  droplet  surface.  Therefore, 
the  change  in  the  predictions  of  droplet  lifetime,  under  the  occurrence  of  thermody¬ 
namic  non-equilibrium  conditions,  would  be  smaller  than  10  %  with  respect  to  the 
thermodynamic-equilibrium  predictions. 

The  saturation  ratio  corresponding  to  the  kinetic  limit  of  superheat  is  reported  for 
the  considered  pressures  in  Fig.  11.  As  the  pressure  is  increased,  the  saturation  ratio 
decreases,  decreasing  thus  the  difference  between  saturation  conditions  and  superheat 
conditions.  This  leads  to  a  ratio  of  enthalpies  of  vaporization  closer  to  unity.  Thus, 
the  influence  of  thermodynamic  non-equilibrium  conditions  drops  to  zero  as  pressure 
increases. 

The  analysis  presented  above  may  be  criticized  because  of  the  somehow  inaccurate 
hypotheses  of  modelization  used.  First,  the  Boltzmann’s  statistics  has  been  used  to 
estimate  the  dynamics  of  the  molecules.  However,  the  Boltzmann’s  statistics  should  be 
accurate  enough  so  that  the  order  of  magnitude  is  not  affected,  and  so  should  be  the 
conclusions  of  the  first  part  of  this  study.  The  molecular  fluxes  are  also  used  to  compute 
the  nucleation  rate.  The  expression  of  the  nucleation  rate  is  a  product  of  the  kinetic 
factor  time  an  exponential  factor.  This  latter  factor  has  the  greatest  influence  on  the 
nucleation  rate  value,  because  as  the  liquid  phase  becomes  more  metastable,  it  increases 
from  a  very  small  value  to  an  almost  infinite  value,®  and  finally  any  inaccuracy  in  the 
estimation  of  the  kinetic  factor  may  be  not  considered  as  dramatic.^’  Second,  our 
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analysis  is  based  on  a  thermodynamic  modelization  of  both  the  bulk  liquid  and  bubble 
gaseous  phase  issued  from  the  Soave-Redlich-Kwong  equation  of  state.  This  cubic 
equation  of  state  has  been  demonstrated  to  lead  to  important  error  in  the  prediction  of 
the  liquid  density  near  the  saturation  line.®  However,  since  first  and  second  derivatives 
of  the  primary  thermodynamic  variables  are  needed  in  the  analysis,  a  relatively  simple 
expression  for  the  equation  of  state  is  needed.  On  one  hand,  one  can  argued  that  this 
may  affect  significantly  the  results,  but  one  must  keep  in  mind,  on  the  other  hand,  that 
the  best  predictions  of  the  droplet  vaporization  times  are  based  on  the  same  kind  of 
thermodynamic  modelizations . 

4  Conclusion 

A  comprehensive  analysis  has  been  conducted  concerning  the  influence  of  thermody¬ 
namic  non-equilibrium  conditions  on  the  droplet  vaporization  process.  The  study  has 
been  restricted  to  the  oxygen/hydrogen  system,  although  the  same  kind  of  analysis 
may  be  conducted  for  other  systems.  First,  the  vaporization  mass  flow  rates  issued 
from  droplet  vaporization  models  based  on  the  classical  thermodynamic  equilibrium 
assumption  have  been  compared  to  the  molecular  fluxes.  It  has  been  shown  that  the 
molecular  dynamics  does  not  constitute  a  limiting  factor  for  the  reaching  of  the  thermo¬ 
dynamic  equilibrium  at  the  interface,  whatever  the  ambient  temperature  and  pressure 
are,  and  as  long  as  the  droplet  diameter  is  above  10  iim.  However,  some  uncertainties 
remain  in  the  expression  of  the  molecular  fluxes,  uncertainties  essentially  contained 
in  the  so-called  accommodation  coefficient.  The  second  part  of  the  analysis  has  been 
devoted  to  the  case  where  molecular  dynamics  is  not  fast  enough  to  sustain  thermo¬ 
dynamic  equilibrium,  i.e.  small  droplet  diameters  or  low  accommodation  coefficients. 
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The  limit  of  occurrence  of  the  resulting  metastable  liquid  phase  have  been  computed, 
both  on  a  thermodynamic  stability  and  on  a  kinetic  limit  viewpoints.  It  has  been 
shown  that,  at  the  kinetic  limit,  the  superheat  liquid  exhibits  properties  such  that  the 
enthalpy  of  vaporization  of  oxygen  is  only  less  than  10  %  different  from  the  saturation 
case.  So  that,  under  thermodynamic  non-equilibrium  conditions,  the  change  in  the 
droplet  vaporization  lifetime  is  expected  to  be  negligible.  In  conclusion,  it  has  been 
demonstrated  that,  as  the  pressure  increases,  first  the  thermodynamic  equilibrium  is 
most  likely  to  be  achieved,  and,  second,  the  influences  of  non-equilibrium  conditions 
become  smaller. 
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Figure  1:  Schematic  diagram  of  the  interface  between  two  bulk  phases. 


Figure  2:  Comparisons  between  equilibrium  net  vaporization  rates  and  molecular  dy¬ 
namics  estimations;  oxygen/hydrogen  system;  p  =  30  atm,  Too  =1000  K,  Do  =  100  pm. 
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Figure  3:  Comparisons  between  equilibrium  net  vaporization  rates  and  molecular 
dynamics  estimations  at  various  pressures;  oxygen/hydrogen  system;  Too  =1000  K, 
Dq  =  100  fim. 


Time,  ms 

Figure  4:  Comparisons  between  equilibrium  net  vaporization  rates  and  molecular 
dynamics  estimations  at  various  ambient  temperatures;  oxygen/hydrogen  system; 
p  =  30  atm,  Do  =  100  pm. 
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Figure  5:  Nucleation  rate  as  a  function  of  the  saturation  ratio  for  various  pressures; 
oxygen/hydrogen  system;  oxygen/hydrogen  system;  T  =  100  K. 
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Figure  6;  Contourplot  of  the  Gibbs  free  energy  in  the  vicinity  of  the  saddle  point; 
Ti  =  110  Pi  =  30  atm,  5  =  5. 


Figure  7:  Enthalpy  of  vaporization  at  the  spinodal  and  kinetic  limit  of  superheat  related 
to  the  enthalpy  of  vaporization  at  the  saturation  conditions;  oxygen/hydrogen  system; 
p  =  20  atm. 


Figure  8:  Enthalpy  of  vaporization  at  the  spinodal  and  kinetic  limit  of  superheat  related 
to  the  enthalpy  of  vaporization  at  the  saturation  conditions;  oxygen/hydrogen  system; 
p  =  30  atm. 
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Interface  Temperature,  K 

Figure  9:  Enthalpy  of  vaporization  at  the  spinodal  and  kinetic  limit  of  superheat 
related  to  the  enthalpy  of  vaporization  at  the  saturation  conditions  as  a  function  of 

tempera ■^•11  TP*  nwcrpn /VivHmorpn  svRt.pm*  n  =  40  aim 


Figure  10:  Enthalpy  of  vaporization  at  the  spinodal  and  kinetic  limit  of  superheat 
related  to  the  enthalpy  of  vaporization  at  the  saturation  conditions  as  a  function  of 
temperature;  oxygen/hydrogen  system;  p  =  50  atm. 
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Figure  11:  Saturation  ratio  as  a  function  of  temperature  at  the  spinodaJ  and  kinetic 
limit  of  saturation;  oxygen/hydrogen  system;  p  =  20,  30,  40,  50  atm. 
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Abstract 

This  paper  presents  a  numerical  simulation  of  liquid-oxygen  droplet  vaporization  and  combustion 
under  high  pressure  conditions  in  hydrogen  quiescent  atmospheres  over  a  wide  range  of  ambient 
conditions,  including  both  subcritical  and  supercritical  regimes.  Theoretical  modelization  has  been 
refined  for  thermodynamic  behavior  and  transport  properties  estimations.  In  the  case  where  gaseous 
water  diffuses  to  oxygen  droplet  surface,  a  simplified  condensation  model  has  been  implemented.  In 
this  paper,  we  focused  our  attention  to  droplet  response  to  ambient  pressure  oscillations  in  terms 
of  vaporization  or  combustion  rate.  Computations  were  carried  out  over  a  wide  range  of  different 
ambient  conditions  for  both  pure-vaporization  and  combustion  cases.  Results  can  be  applied  to 
combustion  instability  studies  of  liquid  rocket  engines. 


Nomenclature 

a  Equation  of  state  constant  in  Eq.(ll) 
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Abstract 

This  paper  presents  a  numerical  simulation  of  liquid-oxygen  droplet  vaporization  and  combus¬ 
tion  under  high  pressure  conditions  in  hydrogen  quiescent  atmospheres  over  a  wide  range  of  ambient 
conditions,  including  both  subcritical  and  supercritical  regimes.  Theoretical  modelization  has  been 
refined  for  thermodynamic  behavior  and  transport  properties  estimations.  Both  pure-vaporization 
and  reactive  cases  have  been  investigated.  In  the  case  where  gaseous  water  diffuses  to  oxygen  droplet 
surface,  a  simplified  condensation  model  has  been  implemented.  In  this  paper,  we  focused  our  atten¬ 
tion  to  droplet  response  to  ambient  pressure  oscillations  in  terms  of  vaporization  or  combustion  rate. 
The  priority  has  been  given  to  the  understanding  of  the  mechanisms  involved  in  the  droplet  response, 
as  well  as  to  the  special  issues  raised  by  high-pressure  environments. 

Nomenclature 

a  Equation  of  state  constant  in  Eq.(ll) 

A  Surface  area 

b  Equation  of  state  constant  in  Eq.(ll) 

c  Speed  of  sound 

Cp  Constant-pressure  specific  heat 

D  Droplet  diameter 

Dij  Binary  mass  diffusion  coefficient 

Dim  Effective  mass  diffusion  coefficient 

e  Specific  total  internal  energy 

Ft  Thermophoretic  force 
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Fy  Viscous  force 

/  Frequency 

hi  Specific  enthalpy  of  species  i 

Kyap  Vaporization  kinetic  coefficient  in  Eq.(21) 

m  Droplet  mass  evaporation  rate 

N  Number  of  species 

p  Pressure 

qe  Energy  diffusion  flux 

qi  Mass  diffusion  flux  of  species  i 

r  Radial  coordinate 

R  Droplet  radius 

Rp  Response  function 

Ru  Universal  gas  constant 

S  Saturation  ratio 

V  Velocity 

Vs  Control  surface  moving  velocity 

V  Total  volume 

Vi  Diffusion  velocity  of  species  i 

Xi  Mole  fraction  of  species  i 

Yi  Mass  fraction  of  species  i 

W  Molecular  weight 

Z  Compressibility  factor 

Greek  Symbols 

a  Thermal  diffusivity 

a{T)  Soave  temperature  function  in  Eq.(ll) 

A  Thermal  conductivity 

Pi  Chemical  potential  of  species  i 

Vi  Stoichiometric  coefficient 

p  Density 

r  Droplet  lifetime 

u)  Acentric  factor 

u)i  Rate  of  production  of  species  i 

Subscripts 

c  Critical  condition 

i  Species  i 

p  Condensed  particle 

r  Reduced  thermodynamic  property 

rf  Reference  fluid 

sat  Saturation  value 

oo  Ambient  condition 


Superscripts 


g  Gaseous  phase 

I  Liquid  phase 

0  Dilute  gas  limit 

Time  derivative 
Fluctuation  quantity 
Averaged  quantity 
*  Dimensionless  quantity 


1  Introduction 

This  paper  is  concerned  with  liquid  oxygen  (LOX)  droplet  vaporization  and  combustion  responses  to 
ambient  pressure  oscillations.  The  ambient  gas  is  composed  of  both  hydrogen  and  water.  In  liquid 
rocket  engines,  propellants  are  injected  into  the  combustion  chamber  as  a  spray  of  droplets  which 
then  undergo  a  sequence  of  vaporization,  mixing  and  combustion  processes.  Each  of  these  processes 
has  been  shown  to  play  an  important  role  in  determining  the  instability  characteristics.  In  this  study, 
special  attention  has  been  given  to  the  droplet  behavior  at  high-pressure  conditions. 

In  the  past  decade,  several  theoretical  works  have  been  developed  to  study  droplet  vaporization 
and  combustion  under  high-pressure  conditions.  Most  of  them  were  numerical  simulations  concern¬ 
ing  vaporization  of  either  hydrocarbon  droplets  in  air,^”^  or  liquid  oxygen  droplets  in  hydrogen 
environments.^”®  By  isolating  the  most  important  processes,  these  works  have  contributed  to  the 
comprehension  of  high-pressure  droplet  vaporization  behavior,  although  some  disagreements  still  ex¬ 
ist  among  predictions,  mainly  due  to  the  transport  and  thermodynamic  modelization  used.  However, 
only  a  few  numerical  simulations  were  devoted  to  oscillating  pressure  conditions.^ 

Droplet  vaporization  and  combustion  have  long  been  identified  as  a  rate-controlling  process  for 
driving  combustion  instabilities  in  propulsion  system.  However,  a  quasi-steady  response  of  the  vapor¬ 
ization  process  to  pressure  oscillations  cannot  sustain  the  combustion  instability.  Many  analyses®”^^ 
have  been  conducted  and  have  shown  that  droplet  vaporization  response  was  small.  Some  of  the  hy¬ 
potheses  used  in  these  models  were  criticized  and  might  have  altered  the  solution.  A  constant  liquid 
temperature  assumption  could  have  kept  the  amplitude  of  vaporization  response  at  low  values.  This 
assumption  is  questionable  when  the  droplet  heating  time  is  longer  than  the  droplet  lifetime,  which  is 
the  case  for  oxygen  droplets  vaporizing  in  hot  hydrogen  atmospheres.  Recently,  Sirignano  et 
have  developed  models  which  showed  that  droplet  response  amplitude  is  sufficiently  large  to  sustain 
the  instability.  In  addition,  the  effect  of  velocity  fluctuation  was  shown  to  be  more  important  than 
the  one  of  pressure  oscillations.  These  latter  models  were  only  applied  to  low-pressure  situations. 

The  present  work  is  concerned  with  the  vaporization  and  combustion  responses  of  liquid  oxygen 
droplets  in  high-pressure  hydrogen/ water  atmospheres.  Since  the  model  is  one-dimensional,  only  the 
effect  of  pressure  oscillations  is  investigated.  The  model  cannot  give  a  definitive  answer  on  the  role  of 
vaporization  in  determining  combustion  instability,  since  the  effect  of  velocity  fluctuation  is  expected 
to  be  significant.  However  it  gives  insights  into  many  important  issues  concerning  the  behavior 
of  the  oxygen/hydrogen  system,  for  both  low-pressure  and  high-pressure  conditions.  In  the  latter 
case,  when  the  vaporization  regime  is  supercritical,  the  droplet  regression,  and  thus  the  vaporization 


mass  flow  rate,  is  classically  computed  through  the  regression  of  a  critical  characteristic  line,  critical 
temperature,  density,  or  composition.  When  interested  in  the  droplet  response  to  ambient  pressure 
oscillations,  this  issue  needs  to  be  addressed  carefully. 

The  numerical  model  is  established  on  a  solid  foundation.  An  important  effort  has  been  devoted  to 
transport  and  thermodynamic  properties  estimations.  A  unified  treatment  of  thermal  conductivity 
and  molecular  diffusion  coefficient  is  presented.  Thermodynamic  properties  are  derived  from  the 
Soave-Redlich-Kwong  equation  of  state  in  a  self-consistent  manner.  In  addition,  a  model  for  treating 
water  condensation  and  its  ensuing  influence  on  droplet  vaporization  behavior  is  developed.^® 

In  this  paper,  we  first  describe  the  theoretical  formulation  with  emphasis  on  estimation  on  trans¬ 
port  and  thermodynamic  properties.  A  condensation  model  is  briefly  discussed,  followed  by  a  sum¬ 
mary  of  the  numerical  methods  used  in  the  present  work.  The  model  is  first  applied  to  vaporization 
of  oxygen  droplets  in  pure  hydrogen  environments  in  both  subcritical  and  supercritical  environments. 
The  reactive  case  is  then  considered. 

2  Theoretical  Formulation 

Figure  1  shows  the  situation  examined  here,  an  isolated  LOX  droplet  vaporizing  in  a  stagnant  hy¬ 
drogen/water  environment  with  superimposed  periodic  pressure  oscillations.  Since  buoyancy  effects 
and  forced  convection  are  ignored,  the  problem  is  governed  by  a  set  of  conservation  equations  in 
spherically  symmetric  coordinate.  The  following  balance  equations  are  written  in  their  most  general 
form.  The  initial  droplet  temperature  is  uniform  and  subcritical  so  that  a  well-defined  interface  in 
thermodynamical  equilibrium  exists  at  time  zero. 

2.1  Conservation  Equations 

Let  us  consider  an  arbitrary  control  volume  Va(t)  delimited  by  a  surface  Aa{t)  moving  in  an  absolute 
referential  at  speed  Va-  The  mass  balance  equation  may  be  written  as  follows: 

^  f  pdV  +  [  p{v-  Va)  ■  dA  =  0  (1) 

dt  yVa(t)  JAait) 

The  momentum  and  energy  conservation  equations  are  simplified  with  the  following  hypothesis:  the 
velocities  involved  are  very  low  so  that  the  kinetic  energy  and  viscous  dissipation  can  be  neglected. 
The  momentum  equation  then  reduces  to  the  following: 

Vp  =  0  (2) 

Thus,  when  the  ambient  pressure  oscillates,  the  droplet  sees  an  environment  with  a  uniform  pressure 
distribution,  but  with  a  time- varying  value.  The  uniform-pressure  hypothesis  is  valid  as  long  as  the 
wave  length  of  the  pressure  oscillation  is  far  greater  than  the  droplet  diameter.  Since  the  wave  length 
is  equal  to  the  ratio  of  the  speed  of  sound  to  the  oscillation  frequency  c//,  the  hypothesis  loses  its 
validity  as  frequency  becomes  very  high. 

The  conservation  of  energy  is  derived  by  applying  the  first  law  of  thermodynamics  to  the  control 
volume 


/  pedV  +  f  peiv  —  u*)  •  dA 
dt  JVa(t)  JAa(t) 


(3) 


The  conservation  of  species  is  written  as 


dt 


/  pYidV+f  pYi{v-Vs)-dA 

JVa(t)  JAait) 


qi  •  dA  +  dji 


(4) 


2.2  Diffusion  Fluxes 

The  diffusion  fluxes  of  energy  and  mass,  qe  and  q{,  in  Eqs.(3)  and  (4)  need  to  be  modeled.  The  Soret 
and  Dufour  effects,  which  account  respectively  for  molecular  diffusion  due  to  temperature  gradient 
and  energy  diffusion  due  to  species-concentration  gradient  are  not  considered  in  the  present  analysis. 
Despite  the  disparity  of  molecular  weight  between  hydrogen  and  oxygen  and  the  high  temperature 
gradients  involved,  these  effects  have  been  proven  to  be  negligible  at  high  pressure  for  oxygen  droplet 
vaporization.^^  For  diffusion  fluxes,  we  derive  the  following  formulas.^® 

9e  =  -Avr-|;?.hi  (5) 

t 

qi  =  pYiVi  =  -pDimVYi  (6) 

Fourier’s  and  Fick’s  laws  are  used  to  approximate  the  thermal  and  molecular  diffusion,  respec¬ 
tively.  The  effective  mass  diffusion  is  obtained  from  the  following  equation  in  terms  of  binary  mass 
diffusivity.^^ 

A„  =  (l-A'i)/f;(J£'i/A,)  (7) 

2.3  Gas-Phase  Combustion  Model 

A  finite-rate  single-step  irreversible  chemical  scheme,  O2  +  2  jff2  — *■  2  H2O  is  considered: 

=  UiWiYo^YH,  (8) 

Combustion  process  is  assumed  to  faster  than  all  the  other  mechanisms  involved.  Nominal  values  of 
A  and  Ea  are  thus  chosen,  so  that  they  do  not  affect  the  solution. 


2.4  Condensation  of  Water 

When  oxygen  reacts  with  hydrogen,  the  main  combustion  product  is  water  vapor.  The  water  vapor 
diffuses  from  the  flame  towards  the  droplet  surface.  Water  is  much  less  volatile  than  oxygen  and 
hydrogen.  Consequently  its  liquid  state  exists  for  higher  temperatures  and  pressures.  Thus  a  non- 
negligible  amount  of  gaseous  water  may  reach  a  low  temperature  region  where  the  water  partial 
pressure  becomes  locally  greater  than  its  saturation  pressure.  Gaseous  water  may  then  experience 
phase  change,  leading  to  formation  of  ice/ water  particles.  In  the  following  analysis  we  consider  that 


gaseous  water  may  experience  only  one  phase  change,  that  is  condensation;  transition  from  liquid  to 
solid  or  direct  transition  from  the  gaseous  to  solid  state  is  not  considered. 

The  condensation  mechanism  involves  an  array  of  complex  processes  such  as  nucleation  and  par¬ 
ticle  growth  by  germination.  In  addition,  particles  may  interact  with  the  gaseous  flow  and  other 
particles.  The  purpose  here  is  not  to  treat  the  entire  problem  induced  by  condensation  because  it 
requires  the  development  of  a  sophisticate  spray  model  which  may  overlap  our  initial  study  set.  All 
the  processes  involved  in  the  condensation  were  analyzed,  which  enables  to  simplify  the  model. 

Since  the  nucleation  process  is  faster  than  any  other  processes  involved  for  realistic  critical  nuclei 
sizes,  condensation  is  assumed  to  occur  instantaneously  when  the  saturation  ratio  {S  =  PH2olPsat) 
is  greater  than  unity.  After  their  formation,  the  condensed  paxticles  interact  with  the  gaseous  flow 
according  to  thermophoretic  and  viscous  forces.  Heavier  particles  are  attracted  to  low  temperature 
regions  due  to  thermophoretic  effect,  based  on  expressions  for  thermophoretic  forces  first  derived  by 
Epstein^^  and  later  corrected  by  Brock. Since  nuclei  are  typically  of  the  order  of  the  molecular 
mean  free  path,  viscous  forces  were  modeled  using  the  Stokes-Cunningham  equation  which  accounts 
for  non-continuum  effects  coincident  with  water  particles  sizes  of  the  molecular  mean  free  path  in  the 
suspending  fluid.  Assuming  that  thermophoretic  and  viscous  forces  axe  the  only  forces  acting  on 
the  particle,  the  equation  of  motion  describing  particle  dynamics  can  be  written  as: 

=  (9) 

Actually  a  terminal  velocity  is  reached  when  the  net  force  acting  on  paxticles  is  equal  to  zero, 
that  is  when  thermophoretic  and  viscous  force  are  equal  and  opposite.  The  time  needed  to  reach  this 
mechanical  equilibrium,  which  expression  may  be  easily  obtained  and  estimated  for  representative 
conditions,  is  by  several  order  of  magnitude  smaller  than  the  other  representative  times,  so  it  is 
possible  to  consider  that  particle  velocity  is  instantaneously  equal  to  its  terminal  velocity  given  by: 

1  dT 

terminal  —  V  0  —  (ifl) 

The  expression  of  0  is  function  of  paxticles  radius,  of  the  ratio  of  thermal  conductivity  for  gaseous 
flow  and  water  particles,  gaseous  viscosity,  and  of  coefficients  used  in  the  expression  of  viscous  forces. 
The  influence  of  thermophoretic  forces  is  rather  small  and  terminal  velocity  is  very  close  to  gas  phase 
velocity.  Furthermore,  since  the  particle  radius  effect  is  entirely  included  in  this  term,  the  overall 
effect  of  condensed  particle  size  is  then  expected  to  be  small  or  even  negligible.  A  full  treatment  of 
water  condensation  process  would  require  a  different  treatment  of  each  condensed  particle  depending 
on  its  size.  Since  each  particle  would  grow  or  decrease  depending  on  the  local  saturation  conditions, 
and  since  each  particle  may  collide  with  other  paxticles,  tracking  each  particle  size  and  position 
would  necessitate  a  sophisticate  modelization.  From  previous  considerations,  such  modelization  is 
unnecessary  and  all  condensed  paxticles  are  assumed  to  have  the  same  radius,  taken  to  be  equal  to 
10“^  pm. 

The  condensation  process  may  then  be  reduced  to  the  addition  of  mass  sink  term,  —mcond,H20 
in  the  balance  equation  of  mass,  Eq.(l),  and  a  source  term,  AHvap,H20  f^cond,H20  in  the  balance 
equation  of  energy,  Eq.(3).  The  clear  parallelism  existing  between  water  condensation  process  and 
oxygen  vaporization  process  leads  us  to  express  these  terms  with  a  method  similar  to  that  used  for 
oxygen  vaporization,  which  is  described  later  in  the  text.  Neglecting  the  contribution  of  Brownian 
motion  of  condensed  paxticles  to  the  pressure,  we  still  can  use  Eq.(ll)  where  p  is  the  gas  density. 


Besides,  we  assume  that  the  total  volume  is  the  sum  of  each  phase  volume.  The  condensed  mass 
volume  is  related  to  its  mass  by  the  particle  density,  which  is  considered  as  independent  of  pressure 
and  temperature. 

2.5  Property  Evaluation 

Transport  and  thermodynamic  properties  play  a  decisive  role  in  determining  droplet  vaporization 
behavior.  The  former  dominates  mass  and  energy  transfer,  whereas  the  latter  affects  the  equilibrium 
conditions  and  energy  needed  for  phase  change.  These  properties  dictate  the  conditions  of  occurrence 
of  supercritical  regimes,  and  as  such,  accurate  predictions  of  thermodynamic  and  transport  properties 
are  needed.  Furthermore,  the  property  evaluation  scheme  must  be  the  same  for  both  liquid  and 
gaseous  phases  to  avoid  consistency  problem  as  the  supercritical  vaporization  regime  is  reached. 

2.5.1  Thermodynamic  Properties 

All  the  thermodynamic  properties  are  derived  from  the  Soave-Redlich-Kwong  equation  of  state.  It 
represents  a  good  compromise  between  computation  complexity  and  physical  accuracy.®  This  equation 
of  state  takes  the  form:^^ 

pRuT  aa 

^  {W  -hp)  W  {W  bp) 

where  a  and  b  are  computed  from  classical  mixing  rules,  a  and  b  accounts  for  attractive  forces  and 
repulsive  forces  between  molecules,  respectively.  All  the  thermodynamic  properties  reqmred  in  the 
present  analysis  can  be  obtained  by  substituting  Eq.(ll)  into  appropriate  thermodynamic  relations.^® 

2.5.2  Transport  Properties 

Transport  properties  determine  the  mass  and  energy  transport  which  dictate  the  interface  conditions 
and  consequently  the  vaporization  rate  and  the  occurrence  of  supercritical  regimes.  The  major 
properties  of  interest  are  thermal  conductivity,  and  mass  diffusivity.  Classical  gas  kinetic  theory 
fails  to  predict  correctly  transport  properties  for  high  pressures.  However  it  is  always  possible  to 
estimate  a  low  pressure  value  and  then  correct  it  by  taking  into  account  high-pressure  effects. 

For  thermal  conductivity  we  use  the  corresponding-state  model  proposed  by  Ely  and  Hanley. 

This  model  is  able  to  predict  transport  properties  over  a  large  domain  of  temperatures  and  pressures, 
from  compressed  liquid  to  dilute  gas.  The  method  requires  the  characteristic  values  of  each  species: 
critical  coordinates  and  acentric  factor.  Details  of  the  procedure  may  be  found  in  the  literature,^^’^'^ 
and  only  a  brief  summary  is  presented  here.  The  mixture  thermal  conductivity  is  divided  in  two 
contributions  as  follows: 

Xmip,T,Yi)  =  Xl{T,Yi)  +  Xl{p,T,Yi)  (12) 

X'^{T,  Yi)  represents  the  internal  contribution  to  thermal  conductivity,  whereas  X'^{p,  T,  represents 
the  translational  and  collisional  parts.  For  a  mixture  the  first  term  may  be  evaluated  by  means  of  an 
empirical  mixing  rule: 

■i  N  N  /i  1 


(13) 


in  conjunction  with  the  Eucken  correlation  modified  for  polyatomic  gases: 

A;.  =  1.32^?(C»,.-j^)  (14) 

where  the  low  pressure  estimations  of  viscosity  and  specific  heat  are  obtained  from  polynomial  ad¬ 
justments  of  existing  experimental  data.^^ 

The  collisional  and  translational  part  \^{p,T,Yi)  is  scaled  to  the  correspondent  property  of  the 
reference  fluid  using  the  extended  corresponding-state  method. 

X'Up,  r,  Y)  =  Kfiprf,  Trf)  TxXy  (15) 

T\  and  X\  are  the  scaling  factor  and  a  correction  factor  respectively.  The  latter  compensates  the 
non-correspondent  behavior  of  the  mixture.  For  the  reference  fluid,  the  coUisional  and  translational 
contribution  of  thermal  conductivity  is  divided  into  three  parts: 

KfiPrf^l'rf)  =  ^rfi^rj)  +  ^Xexc{Prf,Trf) 

+  ^Xcrn{Prfi  Trf)  (16) 

X°f{Trf)  is  the  value  at  the  dilute-gas  limit  for  the  reference  fluid  which  may  be  evaluated  by: 

^  t‘°i 

^Xexc{Prf,Trf)  and  AXcritiprfiTrf)  Correspond  to  the  dense-fluid  correction  and  the  critical  en¬ 
hancement,  respectively.  They  may  be  computed  from  analytical  expressions  fitted  to  experimental 
data.  However,  since  properties  which  diverge  strongly  for  pure  components  diverge  only  weakly  for 
mixtures,^®  and  the  critical  enhancement  of  thermal  conductivity  has  not  been  observed  for  mixtures, 
AXcrit{pTf,Trf)  is  not  implemented  in  the  present  study. 

Our  desire  to  get  a  imified  evaluation  technique  for  mass  diffusivity  for  both  liquid  and  gaseous 
phases  is  inhibited  for  two  reasons.  First,  such  a  technique  does  not  exist  in  the  literature,  second, 
extension  of  existing  modelizations  valid  one  phase  to  the  other  may  lead  to  erroneous  errors.  For 
gaseous  phase,  the  approach  proceeds  in  two  steps.  A  low  pressure  binary  mass  diffusivity  is  first 
obtained  using  the  Chapman-Enskog  theory  in  conjunction  with  the  Lennaxd-Jones  intermolecular 
potential-energy  function.^^  The  coefficient  is  corrected  with  Takahashi  correlation  to  take  into  ac¬ 
count  high  pressure  effects.  Liquid  diffusivity  is  estimated  with  the  Wilke-and-Chang  procedure.  As 
the  supercritical  regime  is  reached,  the  estimation  technique  for  liquids  is  applied  for  temperature 
below  the  mixture  critical  temperature.  These  methods  axe  matched  at  the  interface  by  using  a 
correcting  factor  for  the  liquid  estimation. 


2.6  Interfacial  Boundary  Conditions 


Before  the  surface  temperature  reaches  the  mixture  critical  value,  thermodynamical  phase  equilib¬ 
rium,  chaxacterized  by  a  minimum  of  the  Gibbs  function,  is  assumed  to  prevail  at  the  interface.  If 
surface  tension  is  neglected,  equilibrium  relations  may  be  written  as 

' 

.  =  p'  (17) 

.  i‘’i  =  i4 


The  conservation  of  mass,  species  and  energy  across  the  interface  are 


m  =  p[v-  R)A  =  p{v  -  R)A 

(18) 

m  =  [mYi  =  [rhYi  +  qiAl=ji_ 

(19) 

-  h\) 

(20) 

t 


2.7  Numerical  Methods 

The  set  of  equations,  Eqs  (l)-(4),  is  solved  using  a  fully  implicit  scheme.  Balance  equations  are 
discretized  on  a  time  varying  grid  with  a  finite  volume  formulation.  Two  different  scenarios  occur, 
depending  upon  whether  the  vaporization  regime  is  subcritical  or  supercritical.  In  both  cases  the  grid 
follows  the  time  regression  of  droplet  interface.  In  supercritical  regimes  the  droplet  interface  is  defined 
by  the  critical  characteristic  line  (such  as  temperature,  composition  or  density).  For  each  grid  cell, 
we  compute  the  mass,  the  mass  fraction  of  constituents,  the  internal  energy  and  the  cell  volume.  Cell 
volume  is  derived  from  homothetic  transformation  of  the  initial  grid  and  by  use  of  Eq.(18),  for  tracking 
droplet  regression.  Thermodynamic  quantities  such  as  pressure  and  temperature  are  computed  by 
means  of  a  linearization  technique  using  partial  derivatives  of  temperature,  and  pressure,  respectively 
with  respect  to  the  computed  quantities.  Expressions  of  these  partial  derivatives  are  derived  from 
Eq.(ll).  Details  of  the  numerical  method  are  presented  in  other  publications. 

Let  us  discuss  now  the  special  numerical  procedure  used  in  the  treatment  of  interfacial  conditions 
in  the  subcritical  regime.  The  transfer  of  each  species  involved  in  the  thermodynamical  equilibrium 
between  the  two  phases  is  realized  by  using  the  following  pseudo-phenomenological  or  kinetic-like 
relationship: 

rhi  =  Kyap,i  (21) 

This  expression  is  comparable  to  phenomenological  laws  such  as  Fourier’s  law  for  heat  conduction 
or  Fick’s  law  for  mass  diffusion.  The  flux  is  proportional  to  the  difference  of  potential  via  the 
diffusion  coefficient.  To  ensure  thermodynamical  equilibrium,  the  value  of  the  vaporization  kinetic 
coefficient  is  set  so  that  the  difference  between  the  chemical  potentials  is  as  small  as  needed.  This 
technique  requires  an  implicit  expression  which  is  obtained  by  means  of  a  linearization  technique. 
The  method  does  not  need  any  sub-iterative  process,  but  requires  accurate  expressions  of  chemical 
potential  derivatives.  Once  the  thermodynamical  equilibrium  is  initiated,  the  interface  evolves  on  a 
thermodynamical  equilibrium  trajectory.  This  evolution  is  controlled  by  heat  and  mass  fluxes  at  the 
interface.  As  soon  as  the  interface  reaches  the  mixture  critical  conditions,  this  procedure  is  no  longer 
applied,  and  the  analysis  switches  to  a  single-phase  approach. 

3  Droplet  Vaporization  in  Hydrogen 

In  this  section,  we  consider  the  case  of  pure  vaporization  of  LOX  droplets  in  hydrogen  environments. 
In  this  pure  vaporization  case,  there  is  only  one  length  scale,  that  is  the  initial  droplet  diameter. 
A  dimensional  analysis  shows  that  any  characteristic  time  of  the  process  may  be  related  to  the 


initial  droplet  diameter  by  means  of  a  representative  diffusivity.®  It  follows  that  the  response  of  a 
droplet  to  pressure  oscillations  will  be  identical  for  different  initial  droplet  diameters  as  long  as  the 
product  f  Dq,  or  /i?o,  is  kept  constant.  Figure  2  illustrates  this  point,  the  solid  lines  correspond  to 
a  constant  product  /  Dl  whereas  the  dots  correspond  to  the  computations  which  have  been  carried 
out  for  different  ambient  conditions.  Thus,  a  frequency  of  5,000  Hz  for  a  100  (im  droplet  corresponds 
to  a  frequency  of  555.5  Hz  for  a  300  fim  droplet. 

In  the  present  analysis,  the  initial  droplet  diameter  is  equal  to  100  /fm,  the  amplitude  of  the 
pressure  oscillation  is  5  %.  The  initial  droplet  temperature  is  90  K.  For  this  non-reactive  case  we  are 
interested  in  the  fluctuation  of  the  vaporization  mass  flow  rate  related  to  that  of  pressure,  defined  as 
the  pressure-coupled  vaporization  function  as  follows: 

m  {m  —  rh)frh 
P  ~  [P-  P)/P 

The  Rayleigh  criterion  may  be  simply  stated:  a  wave  will  grow  or  decay  if  sufficient  energy  is  added 
in  phase  with  the  pressure.  If  the  reactive  mechanism  is  assumed  to  occur  without  phase  delay  with 
respect  to  the  vaporization  process,  the  later  criterion  may  be  characterized  by  the  real  part  of  the 
previous  response  function,  which  must  be  positive  to  sustain  instability.  An  estimation  of  this  real 
part  may  be  obtained  by  the  following  formula: 


S(i?,(()l  = 


At+1/2/  •  '  ' 

Jt-1/2/  ^  P 

rt+l/2/  /2  7j 

Jt-imP  it 


(23) 


Let  us  first  consider  the  case  where  the  vaporization  regime  is  subcritical.  In  Fig.  3,  we  present  for 
three  different  pressures  the  time  variations  of  the  droplet  vaporization  mass  flow  rate  fluctuation  m  . 
The  ambient  temperature  is  1500  K,  and  the  oscillation  frequency  is  5,000  Ik.  For  each  pressure,  the 
magnitude  of  the  fluctuation  first  increases,  then  decreases  slowly  until  the  end  of  droplet  lifetime. 
The  amplitude  of  the  fluctuation  increases  with  increasing  pressure.  However,  before  reaching  any 
conclusion,  one  must  look  at  the  phase  angle  between  pressure  and  vaporization-rate  fluctuations. 

In  this  pure  vaporization  case,  the  process  is  entirely  driven  by  transport  phenomena  with  all  the 
characteristic  times  being  related  to  the  square  of  the  instantaneous  droplet  radius.  Furthermore, 
since  the  droplet  response  is  identical  for  a  given  value  of  f  Rq,  it  is  thus  legitimate  to  relate  droplet 
response  to  the  product  f  R^.  This  product  may  be  non-dimensionalized  by  considering  the  liquid 
thermal  inertia  time  a/  =  XiIpCpi.  Figures  4  and  5  show  the  magnitude  and  phase  angle  of  droplet 
vaporization  response  in  terms  of  f  R^fai  and  for  p  =  30  atm.  First,  the  magnitude  and  the  phase 
angle  vary  throughout  the  droplet  lifetime.  However,  each  of  the  presented  curves  collapses  in  the 
later  stage  of  the  droplet  lifetime  on  a  single  curve,  which  is  obviously  a  function  of  the  product 
f  R?  joLi.  As  time  goes,  this  product  decreases  because  of  the  decrease  in  droplet  radius,  so  that  the 
curves  presented  in  Figs.  4  and  5  are  described  from  right  to  left  as  the  droplet  surface  regresses.  In 
the  later  stage  of  the  droplet  lifetime,  the  droplet  response  becomes  a  function  of  the  instantaneous 
product  f  R?  I ai.  This  feature  demonstrates  that  steady  conditions  are  reached  and  the  vaporization 
mass  flow  rate  fluctuates  around  the  steady  conditions  value.  The  amplitude  and  phase  angle  of  the 
response  function  depend  only  on  the  instantaneous  value  of  f  R^fai  for  given  ambient  conditions. 

As  the  magnitude  is  an  increasing  function  of  f  R  jai,  the  phase  angle  decreases,  and  finally 
becomes  less  than  to  —90®,  so  that  the  real  part  of  the  response  function  becomes  negative.  In  Fig.  6, 
we  have  reported  the  real  part  of  the  droplet  vaporization  response,  the  real  part  is  computed  either 


by  means  of  the  two  previous  figures  or  by  use  of  Eq.(23),  both  techniques  gives  consistent  results, 
even  though  the  droplet  response  is  not  perfectly  sinusoidal.  The  real  part  shows  a  maximum  for 
the  product  f  jcci  being  around  300.  The  cut-off  value,  above  which  the  real  part  of  the  response 
function  is  negative,  is  around  3,000. 

A  phenomenological  description  of  the  droplet  vaporization  response  to  ambient  pressure  oscil¬ 
lation  is  difficult  to  obtain  because  many  different  physical  processes  interact  with  each  other  in  a 
fully  coupled  manner.  Furthermore,  a  simplified  formulation  of  the  mechanism,  which  will  simplify 
the  analysis,  does  not  provide  a  satisfying  description.  However,  some  basic  explanations  can  be 
provided  by  considering  the  droplet  surface  conditions  in  terms  of  temperature  and  composition  as 
well  as  the  liquid  thermal  diffusivity  at  the  interface.  First,  although  the  ambient  pressure,  and  thus 
the  ambient  temperature  are  oscillating,  the  interface  temperature  gradient  on  the  gaseous  side  of 
the  interface  is  only  weakly  affected  by  the  pressure  oscillation.  This  is  also  due  to  the  fact  that  the 
characteristic  time  of  the  gaseous  thermal  internia  is  much  smaller  than  the  period  of  the  pressure 
oscillation.  For  similar  reasons,  the  interface  composition  on  the  gaseous  side  is  almost  not  affected 
by  the  pressure  oscillation.  As  the  pressure  oscillates,  the  surface  temperature  oscillates  so  ensuring 
the  interface  thermodynamic  equilibrium.  This  surface-temperature  oscillation  creates  an  oscillation 
of  the  enthalpy  of  vaporization  and,  thus,  of  the  vaporization  mass  flow  rate.  The  average  interface 
temperature  and  pressure  determines  the  magnitude  of  the  oscillation  of  the  enthalpy  of  vaporization. 
The  closer  the  interface  is  to  its  critical  conditions,  the  higher  is  the  oscillation  of  the  enthalpy  of 
vaporization.  This  later  oscillation  controls  the  oscillation  of  the  vaporization  mass  flow  rate.  It  re¬ 
sults  that  the  average  interface  conditions  exerts  a  significant  control  on  the  magnitude  of  the  droplet 
response. 

The  liquid  heat-up  process  influences  both  the  magnitude  and  the  phase  of  the  droplet  vaporization 
response.  The  heat  flux  hitting  the  surface  is  shared  between  the  vaporization  process  and  the  liquid 
heat-up  process.  As  the  liquid  thermal-inertia  characteristic  time,  decreases  with  respect 

to  the  oscillation  period,  the  liquid  phase  absorbs  the  heat  in  its  bulk  volume  and  thus  decreases 
slightly  the  magnitude  of  the  surface  temperature  oscillation.  However,  this  temperature  oscillation 
is  in  phase  with  the  pressure  oscillation.  On  the  contrary,  as  the  product  f  R^jai  increases,  the 
absorption  of  the  heat  by  the  liquid  phase  is  less  easy.  Thus,  only  the  droplet  surface  absorbs  the 
heat  and  with  a  non-negligible  phase  delay. 

The  oxygen/hydrogen  shows  a  very  different  behavior  compared  to  hydrocarbon/air  systems.  For 
n-pentane  droplet  vaporization  in  air,  Hsiao  et  al7  showed  that  the  cut-off  value  of  f  R?  jai  is  much 
smaller,  and  about  200.  The  value  at  which  the  droplet  response  is  maximum  (i.e.,  the  resonant 
value)  is  much  also  smaller.  So  there  is  an  order-of-magnitude  between  both  systems.  Heidmann^® 
predicted  that  the  resonant  frequency,  frequency  at  which  the  maximum  of  the  response  occurs,  for 
hydrocarbon  systems  is  6.5  smaller  than  that  for  oxygen  droplets.  However,  Heidmann’s  work  did  not 
show  any  cut-off  frequency  for  the  oxygen/hydrogen  system  for  the  frequency  range  he  considered, 
probably  because  he  only  considered  values  of  /  R?  jo-i  below  200. 

One  important  difference  between  hydrocarbon  and  oxygen  droplets  is  the  ratio  of  droplet  lifetime 
to  liquid  thermal-inertia  time.  For  hydrocarbon/air  systems,  the  droplet  thermal-inertia  time  is  of 
the  same  order  of  magnitude  as  the  droplet  lifetime.  The  temperature  is  almost  uniform  throughout 
the  droplet  interior  in  the  later  stage  of  the  droplet  lifetime,  whereas  for  the  oxygen/hydrogen  system, 
at  the  end  of  the  droplet  lifetime,  significant  temperature  gradients  still  exist  within  the  droplet.  By 
arbitrarily  increasing  liquid  thermal  conductivity,  the  cut-off  frequency  is  shifted  to  a  smaller  value, 
further  manifesting  this  feature.  This  also  demonstrates  that  simplified  analysis  of  Heidmann  may  not 


be  applied  to  oxygen/hydrogen  system,  since  such  a  theory  postulates  that  liquid  thermal  difFusivity 
is  infinite  and  liquid  temperature  is  uniform.  The  same  observation  was  made  by  Sirignano  et  al.-M 
by  neglecting  the  liquid  interior  temperature  gradients,  the  surface  temperature  oscillations  will  be 
restrained  to  unrealistically  low  values.  In  fact,  the  amplitude  of  oscillation  is  only  slightly  affected, 
but  an  increase  of  liquid  interior  thermal  diffusivity  tends  to  increase  the  phase  lag,  as  mentioned 
previously.  This  causes  the  real  part  of  the  droplet  response  function  to  be  drastically  reduced  and 
eventually  becomes  negative. 

Figure  7,  shows  the  real  parts  of  the  vaporization  response  function  for  three  different  mean 
pressures  in  the  subcritical  regime,  p  =  10,  30  and  50  atm,  respectively.  Two  important  effects  of 
pressure  on  droplet  response  axe  noted.  First,  as  pressure  increases,  both  resonant  and  cut-off  values 
of  f  Oil  are  increased.  Second,  the  maximum  value  of  droplet  response  increases  with  increasing 
pressure.  The  maximum  value  is  even  slightly  above  unity  for  p  =  50  atm.  As  pressure  increases, 
the  droplet  surface  temperature  increases  and  becomes  close  to  its  critical  value.  The  corresponding 
enthalpy  of  vaporization  decreases  and  becomes  more  sensitive  to  temperature  fluctuation.  Conse¬ 
quently,  as  discussed  previously,  the  vaporization  mass  flow  rate  oscillates  in  a  more  important  way. 
In  addition,  for  low  pressure  cases,  the  response  functions  collapse  in  an  earlier  stage  of  the  droplet 
lifetime  with  the  envelope  curve,  showing  that  transient  stage  occupies  a  smaller  paxt  of  the  droplet 
lifetime. 

The  effect  of  ambient  temperature  has  also  been  investigated,  with  the  results  reported  in  Fig.  8. 
As  temperature  increases,  the  maximum  value  of  droplet  response  is  only  weakly  affected;  however, 
the  cut-off  frequency  is  increased.  The  phase  delay  is  reduced,  with  increasing  temperature.  As 
temperature  increases,  the  droplet  lifetime  decreases,  whereas  the  liquid  thermal-inertia  is  almost 
unchanged.  As  mentioned  previously,  this  feature  tends  to  reduce  the  phase  angle  of  the  droplet 
response  function. 

We  now  consider  supercritical  pressures.  An  oxygen  droplet,  initially  with  a  subcritical  tempera¬ 
ture  distribution,  is  injected  into  a  supercritical  environment,  such  that  the  vaporization  regime  may 
transit  from  subcritical  to  supercritical.  The  interface  conditions  evolve  rapidly  so  that  the  liquid 
and  gaseous  phase  become  identical  at  the  interface.  In  previous  works  in  steady  vaporization,  once 
the  supercritical  regime  is  reached,  the  interface  was  defined  based  on  the  locus  of  the  critical  mixing 
temperature.  One  other  definition  based  on  critical  mixing  composition  was  also  used^^  to  account 
for  the  mixing  process.  When  this  latter  definition  is  chosen,  the  critical  interface  tends  to  move 
away  from  the  liquid  core,  mainly  because  of  the  difference  between  mass  and  thermal  diffusion  rates 
(i.e.  Lewis  number  is  not  equal  to  unity).  In  any  case,  it  results  from  the  previous  remaxks  that  the 
droplet  regression  characteristic  time  is  of  the  same  order-of-magnitude  that  the  characteristic  mixing 
time  in  the  gaseous  phase,  either  considering  temperature  diffusion  or  molecular  diffusion,  which  are 
the  two  only  mixing  processes  involved. 

For  low-pressure  condition,  a  large  difference  exists  between  the  physical  characteritics  of  the 
liquid  phase  and  those  of  the  gaseous  phase  issued  from  the  vaporization  process.  The  difference 
in  the  density  enables,  in  spray  modelization,  to  treat  the  vaporizing  droplets  as  source  points  of 
mass  and  sink  points  of  energy.  The  difference  in  the  characteristics  of  both  phases  leads  to  an 
order-of-magnitude  difference  between  the  gaseous  mixing  process  and  the  droplet  regression  due  to 
the  vaporization  process,  so  that  the  gaseous  mixing  process  is  handled  by  the  spray  model.  Thus, 
for  low-pressure  situations,  the  droplet  vaporization  mass  flow  rate  is  an  interesting  and  necessary 
input  to  any  spray  model.  For  supercritical  vaporization  regimes,  no  difference  exists  between  the 
physical  characteristics  of  the  liquid  phase  and  those  of  the  gaseous  phase  issued  from  the  vaporization 


process.  It  results  that  the  characteristic  times  of  the  gaseous  mixing  and  of  the  droplet  regression 
are  of  the  same  order-of-magnitude,  as  mentioned  previously.  Moreover,  the  absence  of  any  disparity 
of  density  between  both  phases  makes  inappropriate  the  treatment  of  the  droplets  as  source  points 
of  mass  in  a  spray  model  under  supercritical  conditions.  This  is  especially  true  when  considering 
non-steady  conditions  of  ambient  pressure  and  velocity,  where  the  droplet  vaporization  response 
is  significantly  aifected  by  the  somehow  arbitrariness  of  the  critical  interface  definition.  However, 
numerical  limitations  prevent  spray  models  to  treat  droplets  without  vaporization  sub-models,  thus 
it  might  be  interesting  to  have  some  qualitative  information  about  the  droplet  vaporization  response 
to  pressure  oscillation  under  supercritical  vaporization  regime. 

In  this  study,  the  critical  interface  is  defined  as  the  critical  density  line.  As  discussed  previously, 
there  is  no  good  definition  of  the  critical  interface.  Different  critical  interface  definitions  lead  to 
significant  change  of  the  phase  delay  between  the  vaporization  response  and  the  pressure  oscillation, 
but  only  slightly  affect  the  magnitude  of  the  response.  The  chosen  definition  leads  to  a  behavior  of 
the  response  which  shares  some  common  features  with  the  sub-critical  vaporization  response,  in  term 
of  phase  angle.  For  supercritical  regimes,  vaporization  mass  flow  rate  is  computed  as  the  mass  flux 
going  through  the  critical  interface. 

Figure  9  shows  the  time  evolution  of  mass  vaporization  rate  fluctuation  for  two  different  ambient 
pressures  of  100  and  200  atm,  respectively.  As  pressure  increases,  the  magnitude  of  the  droplet 
response  increases  also.  For  a  5  %  pressure  oscillation  amplitude,  the  maximum  amplitude  of  mass 
vaporization  rate  fluctuation  is  around  10  %  for  200  atm.  No  visible  effect  on  the  droplet  vaporization 
is  noticeable  as  the  interface  reaches  its  critical  conditions  mainly  because  it  occurs  in  the  very  early 
stage  of  droplet  lifetime.  Figure  10  presents  the  real  part  of  droplet  response  in  terms  of  f  R^/ai  for 
both  pressures.  The  different  curves  correspond  to  various  values  of  /i?o.  Contrary  to  the  sub  critical 
case,  those  curves  do  not  collapse  on  an  envelope  curve.  The  process  becomes  transient  in  nature 
throughout  the  entire  droplet  lifetime.  At  the  beginning  of  the  vaporization  process,  the  response 
function  reaches  a  maximum,  between  1  and  2.  For  p  =  200  aim,  the  phase  angle  becomes  very 
quickly  smaller  than  —90°  so  that  the  real  part  of  the  droplet  response  function  becomes  negative, 
although  the  amplitude  of  the  mass  vaporization  rate  fluctuation  is  higher,  as  shown  in  Fig.  9.  This 
result  is,  however,  very  sensitive  to  the  definition  of  critical  interface  and  must  be  handled  with  care. 

4  Droplet  Combustion  in  Hydrogen 

In  this  section,  the  combustion  process  between  oxygen  and  hydrogen  is  considered.  The  oxygen 
droplet  vaporizes  in  an  atmosphere  composed  initially  of  both  hydrogen  and  water.  Ambient  mass 
fractions  of  both  hydrogen  and  water  are  arbitrary  chosen  to  be  1/2.  For  droplet  burning  in  quiescent 
atmospheres,  the  flame  location  is  driven  by  the  diffusion  of  reactants  outside  the  flame,  and  by  the 
value  of  the  stoichiometric  mixture  ratio,  vo2Wo2Ivh2^H2-  In  the  oxygen/hydrogen  system,  the  high 
values  of  the  diffusion  coefficient  and  stoichiometric  mixture  ratio  render  the  flame  to  be  close  to  the 
oxygen  droplet  surface.  Thus,  temperature  gradients  are  important  at  the  interface,  which  leads  to  a 
very  thin  condensation  zone.  However,  condensation  process  is  not  negligible  since  it  disables  water 
to  reach  oxygen  droplet  surface. 

After  an  initial  transient,  only  oxygen  is  present  at  the  interface:  hydrogen  is  consumed  at  the 
flame,  gaseous  water  condenses  before  reaching  the  droplet  surface  and  condensed  particles  axe  blown 
away  from  the  surface.  Figure  11  demonstrates  this  feature  by  presenting  the  instantaneous  distribu- 


tions  of  species  composition  for  the  three  species  involved.  The  ambient  mean  pressure  is  p  =  30  har. 
The  spikes  in  the  mass  fraction  of  condensed  water  demonstrate  that  the  condensation  region  is 
thin,  and  very  close  to  the  oxygen  droplet  surface.  In  subcritical  gasification  regimes,  the  surface 
temperature  equals  the  oxygen  saturation  temperature  corresponding  to  the  instantaneous  pressure. 
Our  thermodynamic  modelization  is  therefore  not  able  to  handle  computational  cases  where  ambient 
pressure  oscillates  around  the  critical  pressure  of  oxygen.  For  such  cases,  the  gasification  regime  will 
switch  from  a  subcritical  regime  to  a  supercritical  regime  alternately. 

In  this  reactive  case,  it  is  interesting  to  follow  either  gasification  or  burning  response  to  pressure 
oscillations.  The  burning  response  is  defined  in  a  similar  manner  as  the  one  depicted  in  Eq.(22) 
and  with  considering  the  burning  mass  flow  rate,  which  is  the  rate  of  consumption  of  oxygen.  Early 
theories  were  based  on  the  assumption  of  quasi-steadiness  of  the  gaseous  phase,  the  mass  burning 
rate  was  thus  equal  to  the  gasification  rate.  No  vapor  accumulation  between  the  droplet  surface 
and  the  flame  were  taken  into  account.  In  practice,  vapor  accumulation  occurs,  and  consequently 
the  flame  stand-off  ratio  (i.e.,  ratio  of  flame  radius  to  the  instantaneous  droplet  radius)  varies  with 
time.  Especially,  at  the  end  of  droplet  lifetime,  this  ratio  becomes  infinite  during  the  time  needed 
to  burn  the  oxygen  remaining  inside  the  flame.  This  feature  has  also  some  important  consequences 
considering  both  gasification  and  combustion  response  function.  Figure  12  presents  the  time  variations 
of  the  fluctuation  of  gasification  and  combustion  rates  for  three  different  mean  ambient  pressures, 
p  =  10,  30,  40  atm  respectively.  For  the  combustion  case,  the  magnitude  of  the  gasification  response 
function  increases  with  increasing  pressure,  although  it  is  much  smaller  than  for  the  pure  vaporization 
case.  However,  the  main  interesting  characteristic  is  the  behavior  of  the  combustion  response.  As 
one  would  expect,  the  burning  response  follows  the  vaporization  with  a  phase  delay.  This  delay  may 
be  simply  related  to  the  time  needed  for  the  oxygen  vaporized  at  the  droplet  surface  to  reach  the 
flame.  If  the  acoustic  time,  1//,  becomes  much  greater  than  this  characteristic  time,  combustion  and 
gasification  response  should  be  in  phase.  But,  the  most  significant  feature  of  the  combustion  response 
function  is  its  small  amplitude  compared  to  the  gasification  function.  The  diffusion  flame  properties 
may  explain  plausibly  this  point.  As  mentioned  before,  the  location  of  the  flame  is  determined  by 
the  diffusion  of  hydrogen  towards  the  flame.  On  one  hand,  at  the  flame,  the  convection  of  oxygen 
from  the  droplet  surface  by  the  flow  induced  by  the  vaporization  process  is  strongly  influenced  by  the 
pressure  oscillations.  On  the  other  hand,  the  diffusion  of  hydrogen,  which  feeds  the  flame,  is  almost 
not  influenced  by  the  pressure  oscillations.  There  is  thus  an  important  damping  in  the  combustion 
process  provided  by  vapor  accumulation  and  by  the  characteristics  of  the  diffusion  flame.  Nevertheless, 
this  is  a  specific  feature  of  the  one-dimensional  formulation,  the  situation  may  be  completely  different 
under  forced-convection  conditions. 

The  magnitude  of  droplet  combustion  response  is  quite  small,  and  its  phase  angle  is  almost  always 
less  than  —90°,  so  that  the  real  part  of  droplet  combustion  response  does  not  present  any  interesting 
aspects.  Figure  13  reports  the  droplet  gasification  response  for  three  different  subcritical  pressures. 
The  droplet  gasification  response  function  shares  some  common  features  with  the  pure  vaporization 
case.  The  real  part  increases  with  the  average  ambient  pressure.  However,  for  the  range  of  the  product 
/i^  investigated,  neither  maximum  a  of  droplet  response  nor  a  cut-off  value  of  the  factor  fR^fai 
is  observed.  Nevertheless,  for  the  different  cases  simulated,  the  real  part  of  the  droplet  gasification 
response  remains  small. 

For  ambient  pressures  above  the  oxygen  critical  pressure,  the  large  temperature  gradient  induces 
by  the  combustion  process  allows  the  droplet  surface  to  reach  its  critical  conditions  in  the  early 
stage  of  the  droplet  lifetime.  Figure  14  reports  the  time  evolution  of  the  droplet  gasification  and 


combustion  responses  for  two  supercritical  pressure.  The  droplet  response  is  similar  to  the  subcritical 
regime  case.  The  magnitude  of  droplet  combustion  response  is  much  smaller  than  the  corresponding 
gasification  response.  Besides,  the  reaction  process  responds  with  a  more  significant  phase  delay  to 
pressure  oscillations.  One  motivation  of  considering  reactive  cases  was  to  overcome  the  difficulties 
encoimtered  in  the  definition  of  critical  interface  for  supercritical  vaporization  regimes,  by  analyzing 
droplet  combustion  response  instead  of  gasification  response.  The  former  is  insensitive  to  critical 
interface  definition.  However,  since  too  much  difference  exists  between  these  two  droplet  responses, 
comparing  droplet  gasification  response  to  droplet  combustion  response  for  supercritical  regime  does 
not  provide  any  interesting  information,  which  might  have  justified  a  posteriori  the  definition  of 
critical  interface. 

In  Fig.  15,  the  real  part  of  the  droplet  response  is  reported  for  two  different  supercritical  ambient 
pressures.  The  behavior  of  the  droplet  response  is  different  from  the  pure  vaporization  case  at  the  same 
pressure.  Contrary  to  the  non-reactive  case,  after  the  initial  period,  the  droplet  response  becomes  a 
function  of  the  instantaneous  product  fR'^fai,  regardless  of  the  value  of  fR^.  This  may  be  attributed 
to  the  fact  that,  contrary  to  the  pure  vaporization  case,  the  temperature  gradient  at  the  droplet 
surface  is  almost  time-constant,  essentially  because  the  flame  location  varies  only  weakly  during  most 
of  the  droplet  lifetime.  Once  again,  the  droplet  gasification  response  to  pressure  oscillations  appears 
to  be  quite  small. 

5  Concluding  Remarks 

A  comprehensive  analysis  of  LOX  droplet  vaporization  and  burning  response  to  ambient  pressure 
oscillations  in  hydrogen  and  hydrogen/ water  environments  at  both  sub-  and  supercritical  conditions 
has  been  conducted.  Special  attention  was  given  to  evaluation  of  transport  and  thermodynamical 
properties.  A  sub-model  was  developed  to  treat  the  condensation  of  water  vapor  when  it  diffuses 
toward  the  oxygen  droplet  surface.  Both  pure  vaporization  and  combustion  cases  were  investigated. 
The  sensitivity  of  results  to  the  definition  of  critical  interface  has  been  discussed.  In  all  the  cases 
investigated  in  this  study,  the  magnitude  of  the  droplet  vaporization  response  appears  to  be  small 
although  positive.  Since  only  one-dimensional  numerical  simulations  were  carried  out,  only  responses 
to  pressure  oscillation  were  analyzed.  Droplet  response  to  velocity  fluctuations  in  a  convective  envi¬ 
ronment  will  be  treated  elsewhere. 
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Figure  1:  Schematic  of  a  vaporizing  LOX  droplet  in  hydrogen  and  water  environments  and  pressure 
oscillations. 


Fluctuation  of  Mass  Flow  Rate 


Figure  2:  Frequency  function  of  initial  droplet  diameter  for  different  computational  cases. 
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3:  Time  variations  of  fluctuation  of  vaporization  mass  flow  rate  for  different  pressures. 


Figure  4:  Magnitude  of  droplet  vaporization  response  fimction;  p  =  30  atm. 


Figure  7:  Real  parts  of  droplet  vaporization  response  function  for  three  different  mean  pressures; 
subcritical  vaporization  regime. 
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Figure  8:  Real  part  of  droplet  vaporization  response  function  for  three  different  ambient  temperatures; 
p  =  30  atm. 
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Figure  9:  Time  variations  of  fluctuation  of  vaporization  mass  flow  rate  for  different  mean  pressures; 
supercritical  vaporization  regime. 


Figure  10:  Real  part  of  droplet  vaporization  response  function;  supercritical  regime  of  vaporization. 
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Figure  11:  Instantaneous  distributions  of  species  compositions  in  the  entire  field  for  various  times, 
LOXf H2I H-zO  system. 


Figure  12:  Time  variations  of  the  fluctuation  of  gasification  and  combustion  rates  for  different  pres¬ 
sures;  IjOXI H2I H2O  system;  subcritical  gasification  regime. 


Figure  13:  Real  part  of  droplet  gasification  response  function  for  different  pressures;  subcritical 
gasification  regime. 
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Figure  14:  Time  variations  of  fluctuation  of  gasification  and  combustion  rates  for  different  pressures; 
supercritical  gasification  regime. 


Figure  15:  Real  part  of  droplet  gasification  response  function  for  different  pressures;  supercritical 
gasification  regime. 


